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MATHEMATICAL  MODEL  TO  CORRELATE 
FROST  HEAVE  OF  PAVEMENTS  WITH 
LABORATORY  PREDICTIONS 


R.L.  Berg,  G.L.  Guymon  and  T.C.  Johnson 


INTRODUCTION 

Taber  (1930)  stated,  "Freezing  and  thawing  have 
caused  much  damage  to  road  pavements  in  cold  climates, 
but  the  processes  involved  have  not  been  clearly  under¬ 
stood,  and  therefore  some  of  the  preventive  measures 
adopted  have  proved  to  be  of  little  or  no  value.” 

Taber’s  comment  is  still  true  today.  Although  federal, 
state  and  local  governments  have  spent  huge  amounts 
of  money  to  repair  frost  damaged  road  and  airfield  pave¬ 
ments  and  have  made  significant  expenditures  on  the 
development  of  tests  to  determine  the  "frost  suscep¬ 
tibility  "  of  soils,  relatively  little  effort  has  been  expended 
to  develop  a  more  complete  understanding  of  the  physics, 
chemistry  and  mechanics  of  the  frost  heaving  process. 
During  the  development  of  the  mathematical  model  to 
describe  the  frost  heaving  process  that  is  discussed  in 
this  report,  we  learned  that  significant  scientific  advances 
must  still  be  made.  Many  of  these  problem  areas  will 
be  discussed  in  subsequent  sections  of  this  report. 

Within  pavement  systems,  frost  action  is  manifested 
by  heaving  of  the  surface  during  freezing  and/or  a  decrease 
in  the  load  carrying  capacity  of  the  pavement  during 
thawing.  On  roads  severely  affected  by  frost  heave, 
vehicle  running  speeds  must  be  sharply  reduced  to  avoid 
loss  of  vehicle  control  on  the  rough  pavement.  At  air¬ 
fields  which  have  severe  differential  frost  heave  problems, 
it  may  be  necessary  to  close  all  or  a  portion  of  a  runway 
to  aircraft  traffic. 

Conversely,  some  pavements  cannot  support  the 
design  loads  during  "thaw  weakened”  periods,  and  where 
reduced  loads  are  not  imposed  the  deterioration  rates 
of  these  underdesigned  pavements  are  dramatically  in¬ 
creased.  At  some  airfields  it  is  necessary  to  restrict  the 
type  or  weight  of  aircraft  which  may  operate  on  certain 
pavement  features  during  spring  thaw.  At  the  same  time, 
a  number  of  states  impose  load  restrictions  on  some 


highway  pavements  during  the  spring  when  thawing 
occurs. 

The  following  discussion,  contained  in  the  prospectus 
for  this  study,  is  presented  here  to  summarize  the  basis 
for  this  study.  Several  points  in  the  discussion  are  am¬ 
plified  in  subsequent  sections  of  this  report. 

Techniques  to  eliminate  frost  associated  prob¬ 
lems  by  means  of  thermal  insulation  or  soil  mod¬ 
ification  to  prevent  moisture  transfer  have  been 
attempted,  but  are  not  being  used  extensively  be¬ 
cause  they  have  undesirable  side  effects,  are  too 
expensive  or  are  not  universally  effective.  The 
most  widely  applied  technique  for  controlling  the 
detrimental  effects  of  frost  action  is  to  place  non¬ 
frost-susceptible  base  or  subbase  material  to  a 
depth  equal  to  some  fraction  of  the  anticipated 
frost  penetration.  The  availability  of  high  quality 
(non-frost-susceptible)  granular  materials  has 
been  severely  depleted  by  the  construction  of  the 
Interstate  system.  Prudent  use  should  be  made  of 
the  remaining  materials,  conserving  them  for  the 
most  essential  purposes.  To  provide  frost-resis¬ 
tant  pavements  new  sources  of  materials  must  be 
located,  new  standards  of  frost-susceptibility  of 
soils  and  materials  must  be  adopted,  and  tech¬ 
niques  to  beneficiate  unacceptable  materials  de¬ 
veloped. 

Almost  all  frost  heave  susceptibility  criteria 
currently  used  are  based  on  the  particle-size  dis¬ 
tribution  of  the  soil.  Empirical  methods  and  at¬ 
tendant  criteria  based  on  particle-size  distribution 
have  been  widely  used  because  they  are  relatively 
simple  and  rapid.  The  criteria  are  not  consistently 
reliable,  however,  and  soils  meeting  the  established 
criteria  may  be  frost  heave  susceptible,  while 
other  soils  not  meeting  the  criteria  may  actually 
be  non-frost-susceptible.  Criteria  based  on 


particle-size  distribution  are  frequently  misleading 
because  it  is  very  difficult  to  apply  empirical  rules 
to  a  very  complex  phenomenon.  Various  research¬ 
ers  have  identified  a  number  of  parameters,  other 
than  particle-size  distribution,  which  influence 
the  susceptibility  of  soil  and  state  of  stress  in  the 
water;  rate  of  frost  penetration,  overburden  pres¬ 
sure,  permeability,  soil  skeleton  (size,  shape  and 
distribution  of  soil  pores),  cyclic  freezing  and 
thawing,  and  soil  density  are  among  those  factors 
that  have  been  identified  as  influencing  the  severity 
of  frost  heave. 

It  is  unlikely  that  a  consistently  reliable  method 
for  evaluating  or  predicting  a  soil’s  frost  heave 
susceptibility  can  be  based  upon  measurement  of 
only  one  of  the  above  factors,  because  a  number 
of  them  are  interdependent.  A  realistic  laboratory 
assessment  of  a  soil's  frost  heave  potential  requires 
a  test  that  reflects  the  soil’s  behavior  while  it  is 
simultaneously  subjected  to  the  influence  of  all 
these  variables.  To  avoid  the  shortcomings  of  em¬ 
piricism,  several  test  methods  based  on  one  or 
more  of  the  aforementioned  factors,  which  more 
rationally  and  rapidly  evaluate  frost  heave  poten¬ 
tial,  have  recently  been  developed.  The  perme¬ 
ability  test  and  the  heave  stress  test  were  devel¬ 
oped  at  MIT.  The  rapid  frost  heave  test  was 
developed  at  the  University  of  New  Hampshire. 

The  appropriate  freezing  rate  test  developed  by 
Penner  in  Canada  has  also  recently  been  recom¬ 
mended  in  the  literature.  Interpretation  of  early 
results  from  the  first  1  to  2  days  of  freezing  in  the 
CRREL  standard  freezing  test  method  provides 
a  relative  frost  heave  susceptibility  rating.  To 
date  none  of  the  proposed  new  methods  for  rating 
the  frost  heave  potential  of  a  soil  enumerated 
above  has  been  correlated  with  field  performance 
data.  Evaluation  of  the  efficiency  of  these  lab¬ 
oratory  tests  in  categorizing  the  heave  suscep¬ 
tibility  of  particular  soils  requires  that  predicted 
behavior  be  compared  with  observed  behavior  in 
the  field,  where  factors  not  evaluated  in  any  of 
these  laboratory  tests,  such  as  cyclic  freezing 
and  thawing,  may  be  important. 

Research  on  test  equipment  is  essentially  com¬ 
pleted  for  all  the  test  methods  mentioned  above. 
More  laboratory  testing  is  required  to  formulate 
better  concepts  of  the  capabilities  and  limitations 
of  some  of  the  test  methods,  but  correlations  of 
field  performance  with  predictions  by  these  meth¬ 
ods  should  be  undertaken  immediately. 

The  original  purpose  of  the  research  described  in 
this  report  was  to  develop  a  mathematical  model  to  cor¬ 
relate  the  results  obtained  from  the  above  laboratory 
test  methods  with  field  performance.  As  the  study  pro¬ 


gressed  it  appeared  to  the  project  investigators  that  the 
mathematical  model  could  probably  be  used  as  a  design 
tool  to  estimate  the  magnitude  of  frost  heave  for  various 
designs  and  environmental  conditions.  With  further  study 
and  refinement,  the  model  may  also  aid  in  estimating 
the  degree  of  thaw  weakening  and  the  duration  of  thaw- 
weakened  conditions. 

The  two  primary  objectives  of  this  research  were: 

1 .  To  develop  a  mathematical  model  of  frost  heave 
in  terms  of  the  significant  physical,  environmental  and 
climatic  factors,  and  make  a  sensitivity  analysis  of  the 
various  parameters. 

2.  To  develop  a  work  plan,  including  staffing  and 
instrumentation  requirements,  for  correlating  the  results 
of  the  various  laboratory  test  methods  for  assessing  frost 
heave  susceptibility  of  soils  and  the  amount  of  heave 
observed  in  the  field,  making  use  of  the  model  developed 
for  this  purpose. 

The  main  body  of  this  report  addresses  the  first  re¬ 
search  objective,  i.e.  the  development  of  a  mathematical 
model  of  the  frost  heaving  process.  The  second  objective 
was  accomplished  in  late  1976  (Berg  and  Johnson  1976), 
when  an  annual  report  was  submitted  describing  the 
necessary  work.  Important  sections  of  that  report  are 
contained  in  Appendix  A.  Since  that  work  was  accom¬ 
plished,  representatives  of  the  three  agencies  have  dis¬ 
cussed  the  requirements  for  field  and  laboratory  studies 
on  thaw  weakening  of  pavement  systems  in  seasonal 
frost  areas  as  well  as  frost  heaving  of  these  pavements. 

The  FHWA  and  FAA  originally  believed  that  thaw 
weakening  studies  would  be  conducted  subsequent  to 
the  frost  heave  studies.  Because  most  of  the  instrumen¬ 
tation  for  field  observations  on  the  extent  and  duration 
of  thaw  weakening  are  the  same  as  those  for  frost  heave- 
related  studies,  it  was  deemed  desirable  to  conduct  the 
tests  simultaneously.  Appendix  B  contains  a  discussion 
of  the  laboratory  and  field  studies  which  should  be  con¬ 
ducted  to  evaluate  the  thaw  weakening  process. 

ONE-DIMENSIONAL  EQUATIONS  OF  SIMULTANEOUS 
HEAT  AND  MOISTURE  FLUX 

Over  the  last  few  decades  a  large  number  of  reports 
related  to  heat  flux  and  freezing  and  thawing  in  pave¬ 
ment  systems  have  been  published  (Dempsey  1976). 
Johnson  (1952)  lists  more  than  800  publications  on 
these  topics.  Dempsey  (1969),  Christison  (1972),  and 
Berg  (1976)  reference  many  of  the  more  recent  articles 
pertaining  to  freezing  and  thawing  effects  on  pave¬ 
ment  systems. 

Moisture  flux  and  seasonal  changes  in  moisture  con¬ 
tents  in  pavement  systems  have  not  received  as  much 
attention  as  has  heat  flow  in  these  systems.  The  two 
most  prominent  publications  on  this  topic  are  Moisture 
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Equilibrium  and  Moisture  Changes  in  Soils  Beneath 
Covered  Areas  (Aitchison  1965),  the  proceedings  of  a 
conference  held  in  Australia  in  1965,  and  Water  in 
Roads:  Prediction  of  Moisture  Content  of  Road  Sub¬ 
grades,  the  proceedings  of  a  1973  meeting  held  in  France 
by  the  Organisation  for  Economic  Cooperation  and 
Development  (OECD  1973).  Dempsey  (1976)  refers 
to  papers  from  these  two  meetings  as  well  as  several 
other  important  publications. 

Most  literature  on  coupled  heat  and  moisture  flux 
in  soil-water  systems  has  been  published  in  the  last 
decade.  Harlan  (1972)  and  Guymon  and  Luthin  (1974) 
were  among  the  first  to  apply  numerical  models  to 
simultaneous  heat  and  moisture  flux.  An  early  article 
by  johnson  and  Lovell  (1953)  contained  two  figures 
which  summarized  extrinsic  and  intrinsic  factors  af¬ 
fecting  frost  action  in  soils.  Dempsey  (1976)  also  showed 
the  figures  and  we  again  present  them  as  Figures  1  and  2 
of  this  report.  Study  of  these  two  figures  will  illustrate 
the  complex  relationships  which  must  be  considered 
in  developing  mathematical  models  of  frost  heaving. 

In  the  remainder  of  this  section  we  briefly  review  the 
equations  of  heat  and  moisture  flux  in  soil-water  sys¬ 
tems  and  describe  how  the  equations  are  coupled.  Em¬ 
phasis  is  on  freezing  soil  conditions. 

The  theory  and  experimental  verification  of  isother¬ 
mal  moisture  storage  and  movement  in  soils  are  well 
advanced  (Kirkham  and  Powers  1972).  Likewise,  the 
theory  of  soil  thermal  states,  dependent  on  soil  mois¬ 
ture  storage  but  independent  of  soil  moisture  move¬ 
ment,  is  well  developed  (Jumikis  1966,  Luikov  1966). 

It  is  nevertheless  generally  recognized  that  the  hydrologic 
and  thermal  states  of  soil  systems  are  coupled,  par¬ 
ticularly  during  freezing  and  thawing  processes.  In  the 
following  discussion  a  description  of  the  procedure  used 
for  coupling  the  equations  of  heat  and  moisture  flux  is 
discussed  and  a  numerical  analog  for  their  solution  is 
presented. 

With  adequate  computer  capability  it  should  be 
possible  to  realistically  model  the  complex  coupled 
soil  moisture  and  thermal  states  of  soil  systems  involving 
freezing  and  thawing.  The  basis  for  such  a  model  is  our 
current  physically  based  knowledge  of  component  pro¬ 
cesses  combined  with  phenomenological  equations  for 
lumped  processes  not  well  understood,  and  empirically 
based  relationships  for  parameters  arising  in  the  model 
assumptions.  Indeed,  with  the  widespread  availability 
of  computers  there  seems  to  be  little  reason  for  modeling 
only  part  of  the  processes  involved  (e.g.  heat  only)  when 
with  a  little  additional  effort,  more  complete  models  can 
be  developed. 

Philip  and  DeVries  (1957)  presented  a  model  based 
on  the  physics  of  anisothcrmal  moisture  movement  in 


porous  media.  Nielsen  ct  al.  (1972)  reviewed  the  most 
prominent  work  relating  to  anisothcrmal  moisture 
movement.  All  of  this  work  is  primarily  for  soil  water 
systems  where  temperatures  are  above  freezing.  More 
recently,  Harlan  (1972,  1973)  Guymon  and  Luthin 
(1974),  Jame  (1978),  and  Taylor  and  Luthin  (1978) 
have  developed  anisothcrmal  models  of  soil  moisture 
movement  involving  freezing  and  thawing.  The  key 
assumption  involving  the  more  recent  models  is  that 
moisture  moves  in  the  freezing  soil  profile  predominantly 
as  liquid  films  and  obeys  unsaturated  flow  theory. 

Harlan  (1972,  1973)  developed  a  one-dimensional  model 
and  presented  a  finite  difference  numerical  scheme,  as 
did  Taylor  and  Luthin  (1978).  Guymon  and  Luthin 
(1974)  developed  a  finite  element  analog  for  one-dimen¬ 
sional  freezing  and  thawing.  Jame  (1978)  developed  a 
one-dimensional  model  of  a  horizontal  soil  column 
based  on  Harlan’s  work. 

Moisture  Transport 

It  is  well  known  that  Richards’  equation  applies  to 
isothermal  unsalurated  flow  in  ice-free  soil  water  systems 
in  which  the  air  phase  can  be  neglected.  According  to 
Jame  and  Norum  (1972)  and  Harlan  (1973),  Richards’ 
equation  also  applies  to  systems  where  ice  may  be  in 
the  soil  pores  to  some  certain  limiting  condition.  To 
elaborate  upon  this  thesis,  consider  a  single  pore  that 
is  completely  filled  with  water.  As  temperatures  are 
lowered  to  the  ice  nucleation  level  (i.e.  T  <  0°C),  a 
single  crystal  of  ice  commences  to  form  more  or  less  in 
the  center  of  a  pore  because  of  the  distribution  of  free 
energy  in  the  system  (Edlefson  and  Anderson  1943). 

If  the  ice  is  regarded  as  an  air  bubble,  a  situation 
analogous  to  unsaturated  flow  exists.  This  view  permits 
the  use  of  unsaturated  flow  theory  where  the  hydraulic 
conductivity  is  changing  with  respect  to  liquid  moisture 
content.  This  view  will  remain  valid  up  to  some  limiting 
condition  which  will  occur  when  ice  in  soil  pores  has 
developed  to  the  point  that  Darcy’s  law  no  longer 
applies  or  interconnected  unfrozen  water  films  have 
been  disrupted.  This  limiting  condition,  however,  is 
probably  at  a  temperature  much  lower  than  that  usually 
encountered  in  situ. 

Our  proposed  model  is  based  on  the  following  as¬ 
sumptions: 

1 .  Darcy’s  law  applies  to  moisture  movement  in 
both  saturated  and  unsaturated  conditions. 

2.  The  porous  media  are  nondeformablc  as  far 

as  moisture  flux  is  concerned.  It  is  well-known, 
however,  that  as  fine-grained  soils  freeze,  ice 
lenses  form  and  that  the  expansion  of  these 
lenses  together  with  the  expansion  of  freezing 
water  in  pores  causes  deformations  called  heaving. 
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Figure  7.  Extrinsic  factors  influencing  frost  action  (from  Johnson  and  Lovell  1953). 


3.  All  processes  are  single  valued;  i.e.  hysteresis  is 
not  present  in  relationships  such  as  the  soil  water 
characteristic  curve. 

4.  Dissolved  salts  have  a  negligible  effect  on  the 
amount  of  unfrozen  water. 

5.  Water  flux  is  primarily  as  liquid;  i.e.  vapor  flux  is 
negligible. 

Anisothermal  moisture  flux  in  a  one-dimensional 
soil  column  is  governed  by  Darcy’s  law  plus  terms  to 
account  for  vapor  flux  and  thermally  driven  moisture 
and  vapor  (this  flow  law  may  be  termed  the  general 
Darcy  law).  Philip  and  DeVries  (1957),  Nielsen  et  al. 
(1972),  and  Cary  and  Mayland  (1972)  review  these  re¬ 
lationships.  For  purposes  of  this  work  we  assume  that 
the  major  components  of  flux  are  those  due  to  liquid 
driven  by  hydraulic  gradients.  For  moist  soils,  and  it 
is  these  soils  that  are  of  primary  interest,  vapor  transfer 
driven  by  thermal  and  hydraulic  gradients  and  liquid  flux 
driven  by  thermal  gradients  are  generally  several  orders 
of  magnitude  lower  than  liquid  moisture  flux  due  to 
hydraulic  mechanisms.  Since  they  may  contribute  a 
significant  convected  heat  term  and  may  cause  heaving, 
hydraulic  mechanisms  are  the  raison  d'etre  for  including 
the  moisture  state  equation  in  the  first  place. 


A  one-dimensional  moisture  flux  relationship  is 
Darcy’s  law 

p  =  -AfM  (1) 

dx 

where  x  -  the  coordinate  axis,  positive  downward 

v  =  the  velocity  flux  in  the  x  direction,  positive 
downward 

0  =  the  total  hydraulic  head  (0  =  0-x,  where  0 
is  the  pore  pressure  in  length  units,  0  <  0 
for  unsaturated  flow,  and  0  >  0  for  saturated 
flow) 

K  =  7f(0,  T)  the  hydraulic  conductivity 

T  =  temperature. 

Substituting  eq  1  into  the  continuity  equation  applicable 
to  freezing  soils 

-C?,  (2) 

dx  dt 

yields  an  equation  similar  to  Richards’  equation,  i.e. 
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(3) 


where 


d(K  d Mjdx}  _  3/C  _  9£+q 
dx  dx  3 1  1 


where  6  =  9(4/,  T,  t,),  the  volumetric  liquid  water  con¬ 
tent 

Q |  =  the  rate  at  which  liquid  water  is  converted 
to  ice 
t  =  time. 

As  eq  3  now  stands,  it  can  be  solved  numerically  with 
approximate  boundary  conditions;  however,  suitable 
approximations  are  introduced  such  that  the  storage 
and  gravitational  components  of  the  equation  are  cast 
in  terms  of  the  state  variable  4/  which  provides  for  some 
computational  convenience  and  is  required  to  apply  the 
finite  element  method.  It  is  assumed  that  K  and  6  can 
be  approximated  by  relationships  proposed  by  Gardner 
(1958);  i.e. 


K  =  K0  4/  >0 

K  =  KqI(-A^+\)  4><  0 

6  =  60  4>  >  0 

«  =  M^w*3+i)  *<° 


(4) 


(5) 


where  K0  and  C0  are  the  saturated  hydraulic  conduc¬ 
tivity  and  saturated  volumetric  moisture  content 
(porosity),  respectively,  and  Ak  and  Aw  are  parameters 
depending  on  the  soil.  Other  relationships  similar  to 
eq  4  and  5  could  be  used,  e.g.  empirically  determined 
equations  for  particular  soils.  The  Gardner  relationships 
are  used  because  they  were  the  first  such  relationships 
to  be  considered  in  this  research.  They  also  serve  to 
illustrate  the  mathematical  technique  used  in  trans¬ 
forming  eq  3  into  an  equation  containing  a  single  state 
variable.  Since  K  and  6  are  assumed  to  be  continuous 
single-valued  functions  of  i p,  this  in  effect  assumes  that 
the  soil  water  system  has  no  ‘'memory."  Many  soils 
exhibit  hysteresis,  which  is,  in  a  systems  sense,  mem¬ 
ory  of  past  states  of  the  system.  Further,  eq  4  and  5  are 
not  functions  of  the  temperature;  this  will  be  discussed 
following  further  mathematical  development. 

Taking  the  partial  derivative  of  eq  4  with  respect  to* 
and  the  partial  derivative  of  eq  5  with  respect  to  t  yields 

ML=K*M.  4>  <0  (6a) 

dx  dx 

*1=0*^  4/<0  (7a) 

dr  dr 


3AkK04r2 

Mk*3+1)2 

0*= 

(-aw*3+i)2 


(6b) 


(7b) 


Equations  6  and  7  are  substituted  into  eq  3  with  the 
result 


d(K  d4i/dx)  dip  _  q*  W*r>.  (8) 

3*  dx  dt 

where  the  primary  state  variable  is  \p.  The  parameters 
K,  K*,  and  6*  are  functions  of  4/,  and  hence  eq  8  is 
nonlinear.  Parameter  (?|  can  be  specified  a  priori  or  be 
computed  from  suitable  empirical  relationships  to  be 
discussed  subsequently.  (Equation  8  is  strictly  applicable 
to  moisture  flux  in  unsaturated  zones  or  partially  frozen 
soil  and  does  not  apply  to  a  saturated  zone;  to  model 
a  saturated  zone  or  water  table  in  the  one-dimensional 
sense  requires  an  interior  interface  condition  specifying 
the  position  of  the  water  table.)  Before  developing  the 
necessary  coupled  heat  conduction  equation,  appropriate 
auxiliary  conditions  for  eq  8  will  be  presented. 

Two  types  of  boundary  conditions  are  considered: 
specified  boundary  conditions  and  moisture  flux  boun¬ 
dary  conditions.  For  the  upper  boundary, 

4/  -  4/y  x  =  0  t  >  0 

li=Fu  x  =  0  t>0  (9) 

dx 

where  4>\j  is  a  specified  pressure  that  may  be  a  function 
of  time  and  Fy  is  a  known  function  of  time  (e.g.  for  a 
zero  flux  condition,  Fy  =1).  The  lower  boundary  has 
similar  boundary  conditions;  i.e. 


—1 

II 

x  =  L 

t>  0 

i*  =  F, 

dx  L 

x  =  L 

t  >  0. 

(10) 

The  remaining  auxiliary  conditions  are  initial  conditions 
of  the  form 

4'  =  'P0 

0  <  x  <  L 

t  =  0 

(11) 

where  \po  =  4> q(*)  and  is  known  a  priori. 
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Heat  Transport 


A  rigorous  derivation  of  the  heat  transport  equation 
can  be  found  in  Bird  et  al.  (1960).  For  purposes  of  this 
report  a  simpler,  but  correct,  derivation  will  be  presented. 
Most  sources  (e.g.  Meyers  1971)  begin  by  making  the 
deterministic-continuum  assumption,  which  leads  to  a 
partial  differential  equation  with  temperature  as  the 
state  variable  and  with  various  parameters  of  state  such 
as  heat  capacity. 

In  the  model  presented  here,  the  first  concept  em¬ 
ployed  is  that  energy  is  conserved.  Thus,  by  considering 
the  various  rate  processes  involved  in  a  particular  process 
and  by  maintaining  an  energy  balance  on  a  one-dimen¬ 
sional  control  volume,  the  appropriate  heat  equation  is 
obtained.  The  various  rate  processes  that  might  be  con¬ 
sidered  are:  conduction,  convection,  radiation,  heat 
storage,  and  heat  generation  (e.g.  latent  heat  effects). 

In  a  freezing  or  thawing  soil  all  of  the  processes  are 
important.  Moreover,  the  soil  water,  almost  always 
present,  is  a  dilute  solution  containing  dissolved  minerals 
which  affect  the  soil  system’s  thermal  properties.  The 
thermodynamics  of  soil-moisture  systems  is  treated  by 
Edelfson  and  Anderson  (1943)  among  others. 

The  one-dimensional  heat  transport  model  is  based 
upon  the  following  assumptions: 

1 .  The  presence  of  dissolved  salts  does  not  affect 
soil  thermal  properties. 

2.  Freezing  or  thawing  is  an  isothermal  process 
and  latent  heat  generation  can  be  accounted 
for  by  a  “bookkeeping"  process. 

3.  Radiation  heat  transfer  only  occurs  at  the  soil / 
air  interface  and  can  be  accounted  for  in  the 
soil  surface  boundary  conditions. 

4.  Convected  heat  is  primarily  in  the  fluid  form 
and  vapor  convected  heat  is  negligible. 

The  energy  balance  equation  is 

£c<  =  ^t  (12) 

where  £c  is  the  rate  of  heat  conduction  into  the  elemen¬ 
tal  volume,  £v  is  the  net  rate  of  heat  convection  into 
the  elemental  volume,  and  £t  is  the  total  rate  of  heat 
energy  stored  in  the  elemental  volume.  The  individual 
processes  in  eq  12  are  expanded  to  determine  the  net 
heat  conduction  and  convection  from  the  x  direction: 

9|x+Ax  lx 
Ax 

°CW  v  (r-^o)|x+Ax-Cw  v  (r-7'o)|x 

Ax  * 


= fc+fv 


(13) 


where 


q  =  heat  flux, 

Ax  =  the  length  of  the  elemental  volume  element, 
£0  =  a  reference  temperature, 

Cw  =  the  volumetric  heat  capacity  of  fluid 
v  =  the  fluid  flux  in  the  x  direction. 


The  variable  q  may  be  replaced  by  Fourier’s  law: 


(14a) 


where  /CT  is  the  thermal  conductivity  of  the  entire  mass 
of  material.  Differential  equation  14a  is  then  related  to 
the  total  rate  that  heat  is  accumulated  in  the  elemental 
volume,  i.e. 

£t  =  Ca  —  (14b) 

'  a  3 1 


where  Ca  is  the  bulk  volumetric  heat  capacity.  Com¬ 
bining  eqs  13,  14a  and  14b,  dividing  through  by  Ax, 
and  taking  the  limit  as  x  *  0,  results  in 


±  IKj—)  (Cw  vT)  =  C  *L. 
3x  T  3x  3x  a  3r 


(15) 


The  middle  term  in  eq  15,  the  so-called  convection 
term,  may  be  expanded  by  the  chain  rule  and  the  re¬ 
sulting  3v/3x  term  canceled,  because  of  continuity,  to 
yield 


JLik  ALl-C 

3x  Kj  dx  w 


(16) 


The  model  was  originally  developed  using  DeVries 
(1966)  method  to  compute  thermal  conductivity  values 
(Guymon  1976).  Currently  the  program  computes 
thermal  conductivities  using  Kersten’s  (1949)  equations 
for  soils.  Alternatively,  the  thermal  conductivity  for 
each  layer  can  be  specified  in  the  initial  conditions. 
Kersten  (1949)  divided  soils  and  pavement  materials 
into  three  classes:  “fined-grained"  (soils  consisting 
primarily  of  silts  and  clays),  "coarse-grained”  (soils 
consisting  primarily  of  sands  and  gravels)  and  “other” 
(including  flexible  and  rigid  pavement  materials,  in¬ 
sulating  layers  and  soils  with  high  organic  contents). 
Thermal  conductivity  values  for  the  “other”  soil  ma¬ 
terials  must  be  provided  as  input  to  the  program,  and 
the  initial  value  input  for  the  material  is  used  through¬ 
out  the  entire  time  period  of  the  solution.  Kersten’s 
equations,  presented  below,  are  used  for  the  coarse- 
and  fine-grained  soils.  In  the  computer  program, 
thermal  conductivity  values  for  a  partially  frozen  ele¬ 
ment  are  determined  from  the  equations  for  frozen 
soils. 
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For  unfrozen  fine-grained  soils, 

Kt  =  [0,9logw-0.1|  lO001^  (17) 

and  for  frozen  fine-grained  soils, 

Kt  =  0.01  (1O)°-22^d+o.085(10)0OO8'rd  (iv).  (is) 

For  unfrozen  coarse-grained  soils, 

/fT  =  [0.7  log  w+0.4|  10ool7d  (19) 

and  for  frozen  coarse-grained  soils, 

Kj  =0.076(1 O)0-013™ 

+0.032(1 0)00146^' (w).  (20) 

These  are  empirical  equations  where  yd  is  dry  weight 
(lb/ft3 ),  w  is  the  total  moisture  content  (water  plus 
ice  %  dry  weight),  KT  is  thermal  conductivity  (Btu/ft 
h°F). 

Since  the  remainder  of  the  computer  program  uses 
data  in  the  cgs  system  of  units,  dry  unit  weights  are 
converted  from  g/cm3  to  lb/ft3  to  compute  the  thermal 
conductivity,  and  then  the  thermal  conductivity  is 
converted  from  Btu-in./ft2  h  °F  to  cal/cm  h  °C.  Mois¬ 
ture  contents  are  converted  from  volumetric  percentages 
to  percentages  on  a  dry  weight  basis  for  the  thermal  con¬ 
ductivity  calculations. 

A  method  similar  to  that  suggested  by  Kersten  (1949) 
is  used  to  compute  the  volumetric  heat  capacity  of  the 
soil-water-ice  mixture.  The  equation  used  to  calculate 
the  specific  heat  of  the  mixture  is 

cm=  (100cs+wcw+wjCj)/(100+wt)  (21) 

where  cm  =  specific  heat  of  the  mixture 
cs  =  specific  heat  of  the  soil 
w  =  liquid  moisture  content,  %  dry  weight 
cw  =  specific  heat  of  water 
Wj  =  ice  content,  %  dry  weight 
C|  =  specific  heat  of  ice. 

The  volumetric  heat  capacity  Ca  is  determined  simply 
by  multiplying  the  specific  heat  of  the  mixture  by  the 
bulk  density  of  the  soil-water-ice  mixture. 

The  volumetric  latent  heat  of  fusion  is  computed 
from 

L  =*f  Td  Mu)  (22) 


where  L  =  volumetric  latent  heat  of  fusion 
hf  =  heat  of  fusion  of  water 
7d  =  dry  unit  weight  of  soil 
w,  =  total  water  content  (water  plus  ice),  %  dry 
weight 

wu  =  unfrozen  water  content,  %  dry  weight. 

DeVries  (1966  and  1952)  and  others  have  proposed 
additional  empirical  relationships  for  some  of  these 
auxiliary  equations.  Because  they  are  not  essential  to 
the  governing  equation  of  state,  the  numerical  solution 
procedure  developed  here  can  easily  be  modified  to  ac¬ 
commodate  different  auxiliary  equations;  i.e.  they  can 
be  “plugged  into”  the  main  system  model  much  like 
appliances  are  plugged  into  an  electrical  circuit. 

Because  water  storage  and  movement  and  thermal 
states  are  assumed  to  be  coupled  throughout  the  entire 
soil  column,  it  is  necessary  to  specify  boundary  con¬ 
ditions  at  the  same  boundaries  as  used  for  the  moisture 
transport  equation.  This  situation  presents  a  problem 
at  both  ends  of  the  column,  since  at  the  ground  surface 
(i.e.  x  =  0)  heat  is  added  or  lost  by  radiation  and  con¬ 
vection  as  well  as  by  conduction.  Moreover,  the  tem¬ 
perature  or  its  gradient  at  x  =  L  may  not  be  known  at 
the  same  point  that  the  boundary  condition  for  moisture 
flux  is  known.  Both  difficulties  could  be  theoretically 
overcome  by  establishing  two  different  solution  domains; 
however,  such  a  procedure  would  make  programing 
complex.  Accordingly,  these  difficulties  will  be  ignored, 
and  boundary  conditions  of  the  following  type  will  be 
used: 


h? 

II 
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t  >  0 
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(23) 
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t>  0 

*L*Hl 
dx  L 

ii 

* 

t  >  0 

(24) 

where  U  and  L  indicate  upper  and  lower  boundary  con¬ 
ditions,  respectively.  The  upper  and  lower  thermal 
boundary  conditions  must  be  located  at  the  same  posi¬ 
tions  as  the  corresponding  moisture  flux  boundary  con¬ 
ditions.  The  upper  or  lower  specified  temperature  may 
be  a  function  of  time,  and  H  indicates  a  specified  heat 
flux  function  that  may  be  time-dependent.  Additionally, 
cq  1 2  and  certain  of  its  auxiliary  equations  require  ini¬ 
tial  conditions  of  the  following  type: 
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a.  Apparent  heat  capacity. 


Figure  3.  Schematic  representa¬ 
tion  of  two  phase  change  algorithms. 
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and  0  =  (0)o 

0  <  x  <  L 
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© 

(26) 

Phase  Change 

Phase  change  is  the  major  heat  transport  phenomenon 
associated  with  freezing  and  thawing  of  soils.  One  method 
of  including  the  effect  of  phase  change  is  to  use  an 
“apparent  heat  capacity”  over  a  subfreezing  temperature 
range  (Anderson  et  al.  1973,  Johansen  1977).  This  con¬ 
cept  is  shown  schematically  in  Figure  3a.  Since  different 
soils  have  various  curves  in  the  freezing  zone,  Johansen 
(1977)  suggested  using  two  straight  lines  as  shown  in 
the  figure  to  approximate  freezing  within  this  temperature 
range.  Since  moisture  flux  and  formation  of  ice  lenses 
could  cause  the  apparent  heat  capacity  and  computed 
temperatures  in  the  freezing  zone  to  vary  in  an  indeter¬ 
minate  fashion,  a  simpler  procedure  assuming  isother¬ 
mal  freezing  (Fig.  3b)  was  used  in  this  model.  Hromad- 
ka  and  Guymon  (in  prep.)  state  that  the  apparent  heat 
capacity  concept  does  not  converge  near  T  =  0°C, 
and  that  time  steps  required  for  sufficiently  accurate 
solutions  are  very  small. 

We  assume  that  all  of  the  water  available  for  freezing 
does  so  at  a  temperature  somewhat  below  the  freezing 
point  of  pure  bulk  water.  The  freezing  point  depression 
is  an  input  parameter  for  the  computer  program  as  is 


the  amount  of  unfrozen  water.  The  unfrozen  water 
content  depends  primarily  on  the  grain-size  distribution 
of  the  soil  particles  and  can  be  determined  by  any  of 
several  methods  (McGaw  and  Tice  1976,  Anderson  and 
Tice  1972,  Keune  and  Floekstra  1 967). 

A  procedure  similar  to  the  concept  of  “excess  degrees,’ 
(Dusinberre  1961 )  was  used  to  locate  the  freezing  front. 
This  method  requires  a  bookkeeping  process  to  ascer¬ 
tain  the  location  of  the  freezing  front.  It  is  also  similar 
to  the  procedure  used  by  Bafus  and  Guymon  (1976). 

The  amount  of  heat  which  must  be  removed  to  freeze 
an  element  is  determined  from  the  product  of  volumetric 
latent  heat  of  fusion  L  of  the  soil  water  within  an  element 
and  the  volume  of  the  element  V.  The  amount  of  heat 
required  to  change  the  average  temperature  of  an  element 
by  one  degree  is  the  product  of  the  volumetric  heat 
capacity  Ca  and  the  volume  of  the  element.  If  the  first 
product,  i  x  V,  is  divided  by  the  latter,  Ca  x  V,  the 
number  of  “excess  degrees"  E  required  to  freeze  the 
element  is  obtained. 

This  approach  is  illustrated  by  considering  a  common 
point  (node)  A  in  a  discretized  one-dimensional  domain. 
Phase  change  is  not  permitted  to  occur  at  node  A  until 
the  requisite  amount  of  latent  heat  has  been  exhausted 
or  added.  The  prevailing  quantity  of  latent  heat  for 
node  A  is  determined  as  described  above.  The  volume 
of  the  soil,  water  and  ice  mixture  associated  with  each 
point  consists  of  a  finite  element  of  the  mixture  atop 
node  A.  After  each  time  increment  where  a  solution  is 
generated,  the  temperature  of  point  A  is  scanned  to 
determine  whether  or  not  it  has  dropped  below  a  pre¬ 
scribed  temperature  where  phase  transformation  is  as¬ 
sumed  to  occur— the  freezing  point  depression  Td.  If 
the  node  has  not  been  entirely  frozen  and  the  computed 
temperature  has  dropped  below  Td,  then  the  quantity 
of  excess  degrees  evolved  during  the  previously  computed 
time  step  becomes 

f,  =  recomputed  temperature  (27) 

where  £t  represents  the  excess  degrees  removed  from 
node  A  during  the  previous  time  step. 

The  quantity  Et  is  subtracted  from  the  latent  heat  (ex¬ 
cess  degree)  accumulator  for  point  A,  and  if  more  latent 
heat  must  still  be  exhausted,  the  computed  temperature 
at  the  node  is  brought  back  to  Td  prior  to  obtaining  a 
new  solution.  This  is  analogous  to  simulating  phase 
change  as  an  isothermal  process. 

If  the  required  amount,  or  more  than  the  required 
amount  of  latent  heat  has  been  extracted,  node  A  and 
the  volume  of  soil  water  assigned  to  it  are  considered 
to  be  frozen.  The  thermal  properties  of  the  volume 
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common  to  node  A  are  then  updated  based  on  the 
volumetric  proportions  of  soil  water  and  soil  ice  present 
in  the  respective  elements.  If  an  excess  amount  of 
latent  heat  is  removed,  the  residual  is  used  to  calculate 
the  value  of  the  temperature  which  is  now  below  Td. 

Since  the  moisture  content  and  thus  the  volumetric 
latent  heat  of  fusion  and  volumetric  heat  capacity  may 
change  during  each  time  step,  values  of  these  parameters 
are  recomputed  after  each  time  increment.  Other  para¬ 
meters  are  updated  at  a  frequency  specified  in  the  initial 
input  data.  Additional  discussion  of  the  phase  change 
algorithm  is  contained  in  the  section  Evaluation  of  the 
Mathematical  Model. 

The  thawing  process  is  handled  in  a  similar  manner. 
Provisions  are  made  to  determine  which  points  are 
frozen  or  thawed  so  as  to  avoid  refreezing  frozen  areas. 
Boundary  conditions  remain  unaffected  by  the  checking 
routines  that  scan  point  temperatures.  After  phase 
transformation  at  a  point  occurs,  the  latent  heat  ac¬ 
cumulator  assigned  to  that  point  is  recomputed.  Recom¬ 
putation  permits  the  simulation  of  long-term  cyclic 
freezing  and  thawing. 

Coupling  Effects 

The  equation  for  fluid  (water)  flux  in  a  freezing  or 
thawing  soil  adopted  in  this  report  is  eq  8  with  appro¬ 
priate  boundary  conditions  as  discussed  above.  The 
equation  of  heat  flux  in  a  freezing  or  thawing  soil  adopt¬ 
ed  in  this  report  is  eq  16  with  its  associated  boundary 
conditions.  These  two  equations  are  coupled  together, 
especially  in  a  freezing  and  thawing  soil,  because  of  the 
parameters  that  arise  in  the  equations. 

Considering  eq  8  first,  each  term  of  this  equation  is 
coupled  to  eq  16  by  virtue  of  the  parameters  K,  K* , 

0*,  and  Q|,  which  are  dependent  upon  temperature 
while  a  soil  is  freezing  or  thawing.  Similarly,  each  term 
ofeq  16  is  coupled  to  eq  8  because  its  parameters  /?T, 
v,  and  Ca  depend  upon  water  content. 

In  a  highly  dynamic  system,  this  coupling  would 
present  a  formidable  problem.  However,  the  natural 
soil  systems  considered  here  are  highly  damped.  Even 
soil  columns  rapidly  frozen  in  the  laboratory  may  be 
considered  reasonably  damped.  This  is  fortuitous  since 
it  will  allow  the  momentary  decoupling  of  the  governing 
equations  of  state.  Consequently,  one  equation  can  be 
solved  while  holding  its  parameters  constant  for  a  lime 
period  and  then  the  other  equation  can  be  solved  while 
holding  its  parameters  constant.  Parameters  can  be  up¬ 
dated  at  any  arbitrary  frequency  in  the  computation 
procedure  to  provide  the  desired  level  of  precision. 


The  frost  heave  algorithm  consists  of  two  components: 
the  first  is  used  to  compute  the  actual  frost  heave 
during  the  previous  time  increment  and  the  second  com¬ 
ponent  to  define  the  pore  water  pressure  (tension)  at 
the  freezing  front. 

The  first  component  is  comparatively  simple.  The 
porosity  of  each  soil  layer,  and  therefore  for  each  element 
in  the  mathematical  simulation,  is  provided  in  the  input 
data.  At  each  time  increment,  the  heat  available  for 
freezing  additional  water  is  computed.  If  the  volume 
of  water  available  for  freezing  is  greater  than  90%  of 
the  porosity,  frost  heave  occurs  [Dirksen  and  Miller 
(1966)  used  a  value  of  85%  in  their  study] .  The  mag¬ 
nitude  of  this  incremental  amount  of  frost  heaving  is 
determined  by  multiplying  the  amount  of  newly  frozen 
water  in  excess  of  90%  of  the  porosity  by  1 .09.  Thus 
we  assume  that  all  of  the  volume  change  is  in  the  ver¬ 
tical  direction  and  that  the  width  and  depth  of  the 
frozen  increment  are  both  1  cm  in  length.  We  also 
assume  that  all  water  which  enters  the  freezing  element 
travels  to  the  freezing  front;  i.e.  the  moisture  content 
of  the  unfrozen  portions  of  the  element  remains  un¬ 
changed. 

The  second  component  of  the  frost  heave  algorithm 
is  based  on  a  capillary  model  discussed  by  several 
authors,  e.g.  Everett  (1961),  Penner  (1957),  Gold  (1957), 
Williams  (1967)  and  McRoberts  and  Nixon  (1975).  The 
capillary  model  is  based  on  the  premise  that  the  ice/ 
water  interfaces  in  soil  pores  are  curved.  A  pressure 
difference  occurs  at  the  boundary  due  to  the  energy  of 
the  interface.  The  pressure  difference  is 

PrPw  =l£i?L  (28) 

r 

where  P>  =  stress  on  the  ice 

Py,  -  pore  water  pressure 
ojw  =  ice/water  surface  tension  (30.5  dynes/cm) 
r  =  equivalent  pore  radius  of  the  soil. 

The  pressure  on  the  ice  may  be  assumed  equal  to  the 
overburden  stress  and  surcharge  loading  (Williams  1967). 
Thus 

pi=s+it  r>d>  (29) 

/ = i 

where  5  =  surcharge  stress 

7j  =  bulk  density  of  soil  layer  j 
=  thickness  of  soil  layer  /. 
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The  bulk  density  of  a  soil  layer  changes  with  time  as 
water  flows  into  or  out  of  the  layer.  It  also  changes  as 
ice  accumulates  in  the  layer.  In  the  mathematical  model, 
each  element  is  considered  as  an  individual  layer. 

The  equivalent  critical  pore  radius  of  the  soil  is  es¬ 
timated  from  the  soil  moisture  characteristic  curve 
(desorption  portion).  The  empirical  criteria  of  Hoek- 
stra  et  al.  (1965)  are  applied  to  correlate  maximum 
heaving  pressure  developed  by  a  soil  with  its  moisture 
characteristic  curve  using  the  following  reasoning: 

The  maximum  pressure  will  be  developed  in 
order  for  the  ice  interface  to  proliferate  through 
the  smallest  pores.  However,  it  is  clear  that  this 
cannot  be  the  only  criterion;  the  number  of 
small  pores  present  is  another.  If  there  were 
only  a  few  pores  of  a  small  size  present,  the 
freezing  front  could  bypass  them.  The  amount 
of  water  held  in  pores  of  a  particular  size  range 
is  given  by  the  slope. ..[of  the  moisture  character¬ 
istic  curves] . 

They  suggested  that  the  tension  Pc  at  the  point  where 
the  slope  of  the  tension  vs  percent  saturation  curve 
decreased  rapidly  was  an  indicator  of  the  maximum 
heaving  pressure  of  a  soil.  The  same  point  was  used  to 
define  the  effective  critical  pore  radius  in  this  study. 

The  following  equation  was  used  to  determine  the  radius: 


where  rc  -  was  defined  above 

oaw  =  air/water  surface  tension  (72.75  dynes/cm 
at  20°C) 

Pc  =  critical  soil  water  pressure  from  moisture 
characteristic  curve. 

Penner  (1959)  stated  that  the  maximum  overburden 
pressure  to  stop  heaving  appeared  to  be  about  twice  as 
great  as  the  maximum  soil  water  tension  required  to 
stop  heaving.  His  statement  was  based  on  results  from 
laboratory  tests  that  he  had  conducted  and  from  data 
reported  by  Beskow  (1935).  Equation  28  does  not 
suggest  how,  or  if,  the  pressure  should  be  partitioned. 
Based  on  Penner’s  comment,  however,  we  arbitrarily 
doubled  the  effect  of  the  pore  water  pressure  in  the 
equation. 

All  the  referenced  observations  were  included  in  an 
equation  for  computing  the  pore  water  tension  (suction) 
at  the  freezing  front: 

2PW  =  Pr  2-°iZ-  (31) 

rc 

or  by  combining  eq  29,  30,  and  31 


/>w  =  1/2[(S+£  7j  tfj)+0iw  PJ°ivt  1  •  (32) 

/=  1 

The  effects  of  overburden  and  surcharge,  if  any,  are 
taken  into  consideration  by  the  first  term  within  the 
brackets.  This  equation  is  used  to  compute  the  soil 
water  suction  at  the  freezing  front.  Generally,  this 
pressure  “drives”  the  water  movement  in  the  system 
and  it  can  be  considered  a  moving  boundary  condition. 
Once  an  element  is  completely  frozen,  it  is  considered 
that  water  does  not  flow  within  it;  i.e.  no  moisture 
flows  in  the  frozen  zone. 

Using  the  equation  presented  above,  we  may  compute 
soil  water  tension  after  each  time  step.  As  water  flows 
into  the  freezing  element  from  below,  the  degree  of 
saturation  of  the  soil  voids  increases.  When  saturation 
exceeds  90%,  heaving  commences  within  the  element. 

The  amount  of  heaving  is  equivalent  to  9%  of  the  volume 
of  water  up  to  90%  saturation  and  109%  of  the  water 
exceeding  a  degree  of  saturation  of  90%. 

It  is  evident  from  the  above  discussion  that  no  heaving 
occurs  behind  the  freezing  front  in  this  model.  In 
laboratory  studies,  several  authors  have  observed  in¬ 
creased  water  contents  in  the  frozen  soil  (Dirksen  1964, 
Dirksen  and  Miller  1966,  Hoekstra  1966,  Jame  and 
Norum  1976).  Loch  and  Kay  (1978)  reported  that 
laboratory  data  indicate  ice  lenses  form  slightly  behind 
the  freezing  isotherm  in  a  silty  soil.  Taylor  and  Luthin 
(1978)  permit  moisture  flux  in  frozen  soil,  but  the 
laboratory  test  results  that  they  simulated  did  not  ex¬ 
hibit  frost  heave. 

Subsurface  frost  heave  gages  were  recently  installed 
in  a  field  test  section  at  CRREL.  The  soil  was  a  silty 
sand  and  data  from  the  heave  gages  indicate  that  heaving 
continued  at  temperatures  substantially  below  freezing. 
Kinosita  (1975)  also  installed  subsurface  frost  heave 
gages  in  test  sections,  and  his  data  indicate  a  significant 
amount  of  frost  heave  within  the  frozen  zone.  At  this 
time  it  is  not  possible  to  formulate  equations  for  es¬ 
timating  the  amount  of  heave  within  the  frozen  zone. 
Miller*  has  recently  received  funds  for  studying  and 
formulating  his  hypothesis  of  “secondary  heaving” 
and  significant  effort  will  be  expended  by  him  in  the 
next  few  months.  When  a  formulation  for  heaving 
within  a  frozen  soil  is  accomplished  it  may  be  incor¬ 
porated  into  the  existing  model. 

DEVELOPMENT  OF  COMPUTER  MODEL 

The  purpose  of  this  section  of  the  report  is  to  de¬ 
velop  general  solutions  to  the  fluid  flux  and  storage 
equation  and  the  heat  transport  equation. 

*R.D.  Miller,  Cornell  University,  Personal  Communication 
1978. 


So-called  analytical  closed  form  solutions  may  be 
found  for  partial  differential  equations  provided  the 
domain  of  solution,  the  boundary  and  initial  conditions, 
and  the  parameters  are  such  that  a  particular  solution 
can  be  found.  As  a  matter  of  interest  no  such  solution 
has  ever  been  found  for  the  fluid  flux  and  storage  equa¬ 
tion  developed  in  the  preceding  section.  Only  a  limited 
number  of  analytical  solutions  have  been  found  for  the 
heat  transport  equation.  Because  the  domain  of  interest 
is  of  arbitrary  finite  length,  the  boundary  conditions 
variable  in  time  and  space,  and  the  parameters  non¬ 
linear  (i.e.  depending  on  the  moisture  and  temperature 
state),  the  only  possible  solution  is  a  so-called  general 
solution.  General  solutions  imply  both  spatial  and 
temporal  discretization  and  other  approximations  as¬ 
sociated  with  parameters.  Because  the  computational 
effort  is  usually  large  for  most  practical  problems  of 
interest,  the  general  solution  must  usually  be  a  numeri¬ 
cal  solution  using  a  computer.  Although  analog  com¬ 
puters  can  be  used,  most  solutions  are  found  by  digital 
computers. 

Such  an  approach  was  used  in  this  study.  The  model 
was  prepared  for  a  modern  high  speed  computer  using 
FORTRAN  IV  as  a  programing  language.  The  model 
has  been  run  on  a  PDP-1 0  at  the  University  of  Califor¬ 
nia  (Irvine),  on  an  IBM  370/155  in  southern  California, 
and  on  the  Honeywell  system  at  Dartmouth  College. 
Before  the  model  could  be  programed,  however,  the 
partial  differential  equations  had  to  be  reduced  to  a 
numerical  analog  suitable  for  computer  programing. 

The  most  efficient  numerical  analogs  for  solving 
partial  differential  equations  involve  deterministic 
approaches  as  opposed  to  stochastic  methods.  The  two 
basic  deterministic  approaches  are  the  finite  difference 
method  and  the  finite  element  method. 


First,  we  will  consider  the  various  finite  difference 
schemes  investigated.  The  finite  difference  method 
approximates  the  partial  derivatives  found  in  the  boun¬ 
dary  and  initial  value  problem  by  direct  substitution 
of  a  finite  difference  between  solution  end  points  for 
the  continuous  partial  differential  operator.  There  are 
two  basic  approaches  for  doing  this:  the  explicit  dif¬ 
ference  method  and  the  implicit  difference  method. 

The  procedure  will  be  illustrated  by  considering  the 
explicit  technique. 

Let  u  in  eq  33  be  a  single  valued,  finite,  and  con¬ 
tinuous  function  of  x  such  that  u  can  be  expanded  in 
Taylor's  series 

2 

u[x+h)  =  u(x)+hu'(x)+tL-u"  (x)+... 

2 

u(x-h)  =  u(x)-hu'(x)+.~u''  (*)-... 

where  h  =  Ax.  By  assuming  that  fourth  and  higher  powers 
of  h  are  negligible,  the  above  can  be  combined  to  yield 

u(x+h)+u(x-h)  =  2 u(x)+h2u"(x), 
or  rearranging  this  result 

u"(x)=— -0-  =JL  [u{x+h)-2u(x)+u(x-h)\ . 

dx7  x  h 2  (34) 

Since  u  is  also  a  function  of  time  t  (i.e.  u  =u(x,  /)] 
a  suitable  finite  difference  approximation  for  9 u/dt 
must  be  found,  e.g. 


Finite  Difference  vs  Finite  Element  Method 

To  simplify  the  comparison  of  the  two  methods,  the 
partial  differential  equations  were  simplified  to  nor¬ 
malized  symmetrical  parabolic  equations  of  the  type 

9 7u  _  9 u  (33) 

9*2  3f 

so  that  analytical  solutions  could  be  derived  in  order 
to  compaie  with  the  numerical  analogs.  For  the  boun¬ 
dary  and  initial  conditions 

u( 0, t)=a-,u(\,t)=b\t>  0 
u(x,  0)  =  f(x),  0  <  x  <  1 


u(t+k)-u(t)  (35) 

dt  k 

where  k  =  Af. 

Now  if  the  x  space  is  descretized  into  equal  incre¬ 
ments/7  and  the  t  space  is  descretized  into  equal  in¬ 
crements  k,  we  can  adopt  a  standard  subscript  notation 
which  specifies  the  solution  u  at  a  particular  point  / 
in  space  or  a  particular  point/'  in  time;  i.e.  u  -  u (i  j. 
Equations  34  and  35  are  combined  as  a  finite  difference 
analog  for  eq  33  and  the  /',/'  subscript  notation  is  adopted 
to  yield 

U;  |  +  1  =  *_ 

'»  I  O 


(</i+i,r2wi,j+wi-i,i)-"i,j- 


an  analytical  solution  is  easily  obtained  (Meyers  1971 ). 
Such  a  solution  is  used  here  for  comparison  purposes. 


That  is,  u  at  point  /  in  space  and  at  time/'+l  is  ex¬ 
plicitly  determined  by  values  of  u  (at  spatial  points 


/+1 ,  /,  and  /-I )  that  were  computed  in  the  previous 
time  step  j.  For  eq  36  to  converge  (Meyers  1971) 


h 2  2 

Other  finite  difference  schemes  investigated  include  the 
implicit  method 

jK^i-«u)=j[^-("i  +  i.i*i 

+  -i.j  + 1) 

+±(«i+  i.j-^ij+t/i.ij)! 

*2  (37) 

where  /t//t2  < 

These  two  basic  methods  were  solved  using  various 
techniques  consisting  of  initial  boundary  value  problems 
with  derivative  boundary  conditions:  Gauss  elimination, 
Gauss-Seidel  iteration,  and  solution  by  successive  over¬ 
relaxation.  (These  solution  techniques  were  also  com¬ 
pared  to  a  finite  element  analog  for  the  same  equation 
(eq  33)  and  will  be  discussed  in  detail  in  the  next  sec¬ 
tion.] 

The  results  of  this  evaluation  are  qualified  since 
program  efficiency  is  dependent  upon  coding  techniques. 
For  example,  one  programer  may  be  able  to  develop 
a  finite  difference  program  that  is  much  more  efficient 
than  another;  i.e.  efficiency  is  defined  as  computer 
execution  time  vs  solution  precision. 

In  general  the  results  indicated  that  the  finite  dif¬ 
ference  techniques  studied  required  a  substantial  num¬ 
ber  of  iterations  to  reach  a  stable  accurate  solution. 

Time  steps  had  to  be  correspondingly  smaller  than  the 
finite  element  methods  which  did  not  exhibit  stability 
problems.  The  difference  in  computer  execution  time 
was  inconclusive  for  the  one-dimensional  problem 
studied.  Because  of  the  need  for  iteration  in  the  im¬ 
plicit  schemes  and  the  requirement  for  smaller  time 
step  sizes,  it  is  probable  that  two-  and  three-dimensional 
solutions  by  the  finite  element  method  would  be  supe¬ 
rior  to  finite  difference  methods.  The  main  advantage 
of  the  finite  element  method  is  considered  to  be  the 
ability  to  deal  with  complex  and  arbitrary  boundary 
conditions.  Thus,  for  the  one-dimensional  case  it  is 
considered  a  toss-up  as  to  which  method  is  better. 
However,  because  the  model  is  planned  to  be  extended 
to  the  two-dimensional  case,  the  finite  element  method 
is  advantageous  for  the  reasons  mentioned  above. 


Finite  Element  Formulation 

The  heat  equation  and  moisture  flux  equation  have 
the  general  form 

(38) 

where  the  parameters  are  the  heat  equation  para¬ 
meters  or  the  moisture  pressure  state  parameters.  The 
state  variable  u  is  either  the  temperature  or  the  pore 
water  pressure.  The  reason  for  writing  the  coupled 
heat  and  moisture  state  equations  in  the  above  form 
is  to  clearly  indicate  that  the  two  equations  are  identical 
in  form  and  to  facilitate  the  derivation  of  their  finite- 
element  analog. 

We  wish  to  apply  the  Galkerin  method  which  is  a 
special  case  of  the  general  method  of  weighted  residuals. 
To  explain  the  concept,  consider  a  general  operator  A 
which  is  positive  definite.  Suppose  A  operates  on 
the  continuous  function  u  as  follows: 

Au  =  f.  (39) 

Now  we  seek  an  estimate  of  u,  which  we  will  call  b, 
such  that  b approximately  solves  the  above;  i.e. 

Ab^f 

where  b can  be  represented  as  the  product  of  an  ap¬ 
propriate  shape  function  <N>  and  the  value  of  the 
state  variable  at  required  model  points  ;  i.e. 
b  =  <N>  .  We  know  that  the  equation  will  not  be 

exactly  solved  by  bso  that  we  have  a  "residual”  error; 
i.e. 

Ab-f  =  r 

where  r  is  some  small  error.  Since  this  small  unknown 
error  varies  over  the  space  in  which  A  defined,  let  us 
think  of  it  as  a  weighted  average  error  that  we  minimize. 
We  can  minimize  r  as  follows: 

/  </V>T  (Ai/~f)  dx  =  0  (40) 

and  find  b.  The  value  of  b  obtained  will  be  our  best 
estimate  oft;  such  that  r  is  minimum. 

Using  eq  38,  cq  40  becomes 

fL  <A!>t  [*,  *}L+k\dx  =0  (41) 

b  l  dx2  Ox  dt  4J 
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where  L  is  the  total  length  of  the  domain  of  solution 
and  <N>J  is  the  transpose  of  the  weighting  functions 
(shape  functions)  to  be  described  later.  In  eq  41  and 
in  subsequent  equations,  u  is  the  approximation  and 
is  used  to  avoid  the  A  notation.  Befbre  proceeding 
further,  it  is  convenient  to  modify  eq  41  so  that  if  the 
state  variable  u  is  approximated  by  a  linear  shape  func¬ 
tion,  u  =  <N>  \u\  ,  the  first  term  of  the  integrand 

of  eq  41  will  not  vanish.  A  suitable  adjustment  is  made 
by  integrating  the  first  term  of  eq  41  by  parts,  i.e. 

fL  <N>T  k<  $^dx  =k,  fl  <N>J  A  ($iL)dx 

0  1  0X2  1  0  dx\dx) 

=  *.  (<N>r  j"  L-/l  *L  d.<.N>\dx), 
\  dx  0  0  dx  dx  ' 

(42) 


The  first  term  in  the  result  is  associated  with  the  boun¬ 
dary  conditions  where  du/dx  =  0  on  the  boundary  if 
u  is  specified.  If  u  is  not  specified,  du/dx  -*  0;  hence 
it  is  called  a  natural  boundary  condition  since  it  is 
automatically  satisfied.  Ignoring  the  boundary  con¬ 
ditions  for  the  time  being  and  substituting  the  second 
term  of  eq  42  into  41  and  multiplying  through  by  -1 . 
we  have 


M 


£ 
m  =  1 


a<y.>L  <n>7  a. 

dx  dx  dx 


+k3  <N>J  —~k4  <N>J  dx  =  0  (43) 


where  the  solution  domain,  0  <  x  <  L,  has  been  divided 
into  M  finite  elements  of  length  ?.  Substituting  the 
shape  function  u  =  <N>  into  eq  43  yields 


M 

£ 


m  ~  1 


d  <N>J  9  <N>  (y) 
dx  3x  '  » 


+k2<N>T  *-<»>  ft/} 

2  dx  1  ’ 

+/?3  <N>J  <N>  -k4  </V>T)  dx  =0 

(44) 

where  the  dot  in  {«}  indicates  differentiation  with 
respect  to  time  t. 

In  the  finite-element  method,  a  solution  of  the  prob¬ 
lem  is  obtained  by  specifying  an  approximation  of  the 
state  variable  in  each  finite  element  of  the  domain  such 
that  1 )  the  state  variable  is  continuous  throughout  the 


entire  domain,  2)  the  first  derivative  with  respect  to 
space  exists  in  each  element  (but  not  necessarily  on  the 
element  interfaces),  and  3)  the  boundary  conditions 
are  satisfied.  These  last  three  conditions  are  called 
“admissibility  requirements.”  First,  consider  the  solu¬ 
tion  domain  divided  into  M  finite  elements  as  follows: 

Elements 

m  =  1  m  =  2  m- 3  m- 4 

0 - -0 - 0 - 0 - 0 - * 

x  =  0  x 

1  3  5  7  9 

Nodes 

Consider  the  mth  element  as  follows: 


Interface 

Auxiliary 

Interface 

Node 

Node 

Node 

2m  -  1 

2m 

2m  +  1 

u 

0 

\) 

x2m  -  1 

1 

x2m 

*2/77  +  1 

£ 


where  x  is  the  global  coordinate  and  £  is  the  local  coor¬ 
dinate.  An  auxiliary  node  is  included  and  will  be  ex¬ 
plained  shortly.  The  relationship  between  coordinates 
is  given  by 

The  length  of  the  mth  element  is 
®  ~x7m  +  l'*2m  -  1 


and  node  2m  is  at  the  midpoint  of  the  element.  We  wish 
to  approximate  the  state  variable  u  by  a  second  order 
polynomial  (quadratic  shape  function)  as  follows: 


w  =  s+6£+c£2=  |1  ££2) 


a 

b 

c 


(45) 


The  unknown  coefficients  {a,  b,  and  c)  can  be  deter¬ 
mined  if  there  are  three  independent,  simultaneous 
equations  to  solve;  hence,  the  purpose  of  the  midpoint 
node: 


u2m-\  _l 1  0  0] 
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u2m  =  [1  e/2  C2/4l  b 

C 

a 

^m  +  I  =  |1K82]  b 
c 

which  can  be  formed  into  a  3x  3  matrix 
u2m  - 1  10  0  a 

u2m  =  1  8/2  £2/4  b  . 

u2m  +  1  1  8  82  C 

Taking  the  inverse  of  the  3x3, 

a  10  0  u2m  -1 

b  =  -3/8  4/8  -1/8  u2m 

c  2/e2  -4/e2  2/e2  w2m+1 

and  substituting  this  into  eq  47  gives 

«=  [1*«21  1  0  0  «2m_, 

-3/8  4/8  -1/8  tf2w  . 

2/82  -4/82  2/B2  u2m  +  1 

(46) 

By  matrix  multiplication  the  above  can  be  written  as 


u  =  l|e2-38?+2|2,  4e?-4|2,  -e*+2£2l  u2m_i 

e2  "2  m 

u2m  +  1 

(47) 

which  can  also  be  written  as 

u  =  <N>  jw}  (48) 

where  u  is  the  state  variable  in  terms  of  the  nodal  dis¬ 
placements  |i/|  ,  <N>  is  a  1  x  3  matrix  which  is  a  func¬ 
tion  of  8  only,  aid  |iz|  is  a  column  matrix  of  the 
nodal  state  variable  values. 

Returning  to  eq  44,  it  is  integrated  term  by  term  to 
obtain  the  desired  result  which  consists  of  a  matrix 
equation  of  the  form 

U!  H  +1  p]  u  =  \f)  (49) 

which  is  derived  in  detail  in  Appendix  C.  This  element 
system  equation  is  appropriately  assembled  with  all 
other  system  equations,  as  is  indicated  in  the  appendix, 
to  form  the  system  equation 

(si  H  M  =  H  (50> 

where  [S]  is  a  square  banded  nunsymmetrical  matrix, 
ju  |  is  a  column  matrix  of  nodal  unknowns,  ju}  is 
the  ordinary  time  derivative,  and  is  a  column 
matrix.  The  )5)  matrix  is  determined  by  the  geometry 
of  the  discretization  in  space  and  the  conductivity 
parameter.  [/’]  is  determined  by  spatial  discretization 
and  the  capacitance  parameter.  [  is  determined 
by  the  ice  sink  term  and  the  boundary  conditions.  As 
is  indicated  in  Appendix  C,  special  computer  storage 
of  eq  50  is  possible  to  minimize  storage  requirements. 
Equation  50  represents  a  system  of  ordinary  equations 
which  must  be  solved  with  respect  to  time. 

Time  Domain  Solution 

Once  system  matrices  have  been  constructed  in  the 
finite-element  formulation  of  a  parabolic  partial  dif¬ 
ferential  equation,  the  problem  becomes  one  of  solving 
the  linear  first  order  differential  equations  given  by  eq 
50.  The  Crank-Nicholson  method  of  solution  is  util¬ 
ized  here  and  consists  of  an  arithmetic  mean  of  the 
derivative  of  the  state  variable  at  the  beginning  and  end 
of  each  time  step  to  move  the  solution  ahead  in  time.  Thus 

Mv+1  =  Mv+f  (5d 

where  v  is  the  time  step. 


15 


Premultiplying  by  the  [P\  matrix,  this  becomes 

[p\  {«}v^ 

([ P\  MW*!  {^v+’)-  (52) 

But 

[P\  jw|  =  jF|-[Sj  |t/}.  (53) 

Therefore, 

[p\  \u\"'  =  [p\  Hv+f 

(|F|-|S|  |»|>.|F(-|SI  H”') 

[S|  +  f[/>]  |</}v+1 

A t 

■(i-iw)  i"r*2ifi 

Since  |u|  v  would  be  the  initial  conditions  of  a  prob¬ 
lem  or  previously  computed  values  of  the  state  variable, 
the  above  equation  can  be  reduced  to 

|H'|  Mv+1  =1*1  Mv+2  \f] 

[W\  |«{  V+I  =  {g(  (54) 

where  V+I  is  computed  by  Gaussian  elimination. 

EVALUATION  OF  THE  MATHEMATICAL  MODEL 

Evaluation  of  the  model  is  being  accomplished  in 
two  phases:  1 )  study  of  the  accuracy  of  the  finite 
element  approach  as  applied  to  problems  where  the 
governing  equations  of  state  are  known  and  exact  solu¬ 
tions  are  available,  and  2)  study  of  the  validity  of  assump¬ 
tions  necessary  to  extend  the  model  to  problems  where 
the  equations  of  state  and  their  solutions  are  not  known. 
Evaluation  of  the  model  has  proceeded  step  by  step. 

In  our  work,  we  have  verified  by  demonstration 
rather  than  by  theoretical  analysis.  Initial  and  boun¬ 
dary  conditions  for  which  analytical  solutions  are 
available  have  been  used  to  evaluate  and  verify  portions 
of  the  model.  Comparisons  with  laboratory  experiments 
have  been  made  to  evaluate  and  refine  other  portions 
of  the  model.  Initially  heat  and  moisture  flux  portions 
of  the  model  were  “decoupled”  to  evaluate  each  por¬ 
tion  individually.  Decoupling  was  accomplished  by 


choosing  heat  or  moisture  flux  parameters  which  ef¬ 
fectively  caused  no  moisture  flux  when  evaluating  the 
heat  flow  portions  and  vice-versa.  The  following  sections 
initially  discuss  the  decoupled  evaluations  and  lastly  the 
coupled  model  using  both  laboratory  and  field  test 
results  for  comparison. 

Heat  Flux 

The  heat  flux  portion  of  the  model  was  compared 
with  analytical  solutions  for  two  different  problems  in 
which  phase  change  does  not  occur.  In  both  cases  the 
convective  heat  transfer  contributions  have  been  set 
to  zero,  and  the  moisture  flow  part  of  the  coupled 
equations  has,  in  effect,  been  neglected.  In  these  com¬ 
parisons,  the  volumetric  heat  capacity  was  0.582  cal/cm3 
°C  and  the  thermal  conductivity  was  16.05  cal/cm  h 
°C,  typical  values  for  a  silty  soil.  First  the  model  was 
compared  to  an  analytical  solution  for  the  transient 
temperature  distribution  in  a  homogeneous  semiinfinite 
body  initially  at  a  uniform  temperature  with  a  step 
change  of  20°C  applied  to  the  upper  surface  of  the 
body.  Figures  4  and  5  show  the  maximum  deviation 
of  the  finite  element  method  (FEM)  solution  (using 
linear  shape  functions)  from  the  analytical  solution  for 
several  different  combinations  of  element  lengths  and 
time  increments  in  elapsed  times  of  1  and  10  hours.  In 
Figure  4  the  time  increment  is  0.1  h  while  the  element 
length  (uniform  over  the  entire  spatial  domain)  is  varied. 

A  progression  toward  smaller  elements  initially  decreases 
the  maximum  error,  but  after  some  apparently  optimum 
element  length,  approximately  2  cm,  the  trend  is  reversed 
and  larger  errors  again  occur.  It  is  important  to  note, 
however,  that  these  errors  occur  at  nodes  just  below 
the  forced  upper  boundary  (column  4  in  Table  1)  and 
that  the  error  decreases  quickly  with  depth.  Table  1 
illustrates  that  the  accuracy  at  an  arbitrarily  chosen 
depth  of  10  cm  steadily  improved  as  the  element  length 
decreased,  although  the  maximum  error  increased  be¬ 
tween  element  lengths  of  2.0  to  1 .0  and  1 .0  to  0.5  cm. 

As  indicated  in  Figure  5,  shrinking  the  time  increment 
decreased  the  maximum  error,  although  the  rate  of  im¬ 
provement  became  insignificant  rather  quickly. 

Table  2  provides  another  perspective  regarding 
errors  associated  with  time  discretization.  The  fourth 
column  shows  depths  at  which  the  maximum  errors 
occurred.  This  depth  increased  as  the  maximum  error 
decreased,  indicating  that  the  initial  error  was  propagated 
downward  through  the  profile  in  successive  time  increments. 
The  fifth  column  contains  errors  at  a  depth  of  1 0  cm. 

The  relatively  small  error  corresponding  to  a  time  in¬ 
crement  of  0.5  hours  is  deceptive.  Since  a  large  negative 
error  exists  at  a  depth  of  8.0  cm  and  a  large  positive 
error  exists  at  1 2.0  cm  (these  data  are  not  shown  in 
this  report),  the  error  observed  at  10.0  cm  falls,  by 
coincidence,  near  the  intersection  of  the  FEM  and 
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Figure  4.  Maximum  errors  after  1  and  10  hours  of  simulation 
for  various  element  lengths  and  a  time  increment  of  0. 1  hr. 


Figure  5.  Maximum  errors  after  1  and  10  hours  of  simulation  for  various 
time  increments  and  element  lengths  of  2  cm. 


Table  1.  Errors  after  one  hour  for  various  element  lengths 
and  a  constant  time  increment. 


Element 

length 

(cm) 

Time 

Increment 

(h) 

Maximum  error * 
in  profile 
(°C) 

Depth  at 
which  max. 

error  occurs 
(cm) 

Error  at 

10.0<m 

depth 

(°C) 

10.0 

0.1 

-1.952 

10.0 

-1.952 

5.0 

0.1 

-0.638 

5.0 

-0.537 

2.0 

0.1 

-0.126 

6.0 

-0.084 

1.0 

0.1 

-0.346 

2.0 

-0.021 

0.5 

0.1 

-2.448 

0.5 

-0.005 

♦Error  =  analytical  solution  •  FEM  solution. 
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Table  2.  Errors  after  one  hour  for  various  time  increments 
and  a  constant  element  length. 


Depth  at 


Element 

rime 

Maximum  error ’ 

which  max. 

Error  at 

length 

increment 

in  profile 

error  occurs 

10.0  cm 

(cm) 

(h) 

ro 

(cm) 

rc) 

2.0 

1.0 

-8.610 

2.0 

-0.813 

2.0 

0,5 

+  3.400 

2.0 

+0.088 

2.0 

0.25 

-0.807 

4.0 

-0.131 

2.0 

0.1 

-0.126 

6.0 

-0.084 

2.0 

0.05 

-0.122 

6.0 

-0.074 

2.0 

0.025 

-0.120 

6.0 

-0.084 

2.0 

0.01 

-0.120 

6.0 

-0.084 

•Error  =  analytical  solution  *  FEM  solution 


analytical  solutions.  This  simple  error  analysis  provides 
some  insight  into  the  performance  of  the  model. 

However,  a  root  mean  square  error  analysis  on  the  ac¬ 
curacy  of  the  finite  element  method  with  second- 
order  shape  functions  has  been  performed  for  a  similar 
problem  and  provides  a  different  perspective  of  com¬ 
putational  errors.  In  the  finite  element  solutions,  the 
overall  accuracy  improved  rapidly  as  the  spatial  in¬ 
crement  was  made  smaller,  until  the  root  mean  square 
error  was  essentially  zero.  A  similar  trend  was  observed 
when  the  size  of  the  time  increment  was  decreased, 
although  in  both  cases,  the  error  near  the  forced  boundary 
condition  may  have  been  considerable.  Data  also  indicate 
that  second  order  shape  functions  use  less  computer 
run  time  in  achieving  a  given  level  of  accuracy  than  do 
linear  shape  functions. 

It  is  important  to  note  that  a  surface  temperature 
step  change  of  20°C  is  a  severe  test  of  the  model,  and 
that  a  step  change  of  this  magnitude  will  rarely  be 
observed  in  actual  road  or  airfield  pavements.  None¬ 
theless,  the  attainable  accuracy  is  well  within  engineering 
requirements.  In  addition,  our  studies  have  shown 
that,  for  all  combinations  of  spatial  and  temporal 
discretization,  a  smaller  surface  temperature  step  change 
caused  correspondingly  smaller  errors;  in  fact,  if  the 
magnitude  of  the  step  change  were  decreased  by 
one-half,  the  new  error  would  be  significantly  less  than 
one-half  the  previous  error. 

The  heat  flow  portion  of  the  finite  element  model 
was  also  compared  with  the  analytical  solution  to  the 
problem  where  the  surface  temperature  varied  sinusoidally. 
The  temperature  of  the  semi-infinite  body  was  initially 
uniform  at  the  mean  temperature  of  the  sinusoid. 

Again,  phase  change  was  not  considered.  The  sinusoid 
was  approximated  in  the  computer  simulation  by  a 
series  of  step  changes,  occurring  at  each  discrete  time 
interval.  The  overall  error  decreased  with  time  for 
both  the  linear  and  second-order  shape  functions,  al¬ 
though  the  error  at  the  node  directly  below  the  surface 


occasionally  increased  slightly.  The  largest  errors  oc¬ 
curred  after  the  largest  step  changes;  with  constant 
time  increments  this  happened  whenever  the  argument 
of  the  sine  function  was  a  multiple  of  n.  The  overall 
accuracy  of  the  solution  was  especially  sensitive  to  the 
size  of  the  time  increment  because  1 )  smaller  time  in¬ 
crements  cause  smaller  step  changes  and  better  represen¬ 
tation  of  the  sinusoid,  and  2)  the  Crank-Nicholson 
method  leads  to  better  convergence.  Again,  second- 
order  shape  functions  appeared  to  be  more  efficient 
(required  less  computer  run  time  for  a  given  level  of 
accuracy)  than  linear  shape  functions. 

The  phase  change  algorithm  is  an  extension  of  the 
heat  flow  problem;  again  the  convective  contribution 
was  effectively  set  to  zero  for  this  evaluation.  Table 
3  contains  the  differences  between  solutions  obtained 
from  the  modified  Berggren  equation  and  those  ob¬ 
tained  using  the  FEM  model  (with  linear  shape  func¬ 
tions)  for  seven  representative  combinations  of  time 
increment  size  and  node  spacing.  Thermal  properties 
of  the  soil  used  in  this  evaluation  are  slightly  different 
from  those  used  previously  and  are  shown  in  Table  3. 
The  ratio  of  latent  heat  accumulated  to  the  total  heat 
extraction  necessary  for  complete  freezing  was  used 
to  approximate  the  position  of  the  freezing  front 
within  the  element.  Data  indicate  that  decreasing  the 
size  of  the  element  does  not  necessarily  improve  the 
accuracy  of  the  solution.  There  are  two  reasons  for 
this:  1 )  the  algorithm  is  not  adequate  if  the  rate  of 
frost  penetration  is  such  that  the  freezing  front  ad¬ 
vances  more  than  one  element  in  a  single  time  incre¬ 
ment,  and  2)  the  first  part  of  the  two-step  algorithm 
may  yield  physically  unrealistic  temperatures  if  the 
time  increment  is  too  large  relative  to  the  element 
size.  A  further  complication  is  that  the  accuracy  of 
the  solution  docs  not  necessarily  improve  if  the  time 
increment  size  is  decreased  while  the  element  length 
is  held  constant.  The  reasons  for  this  behavior  are  not 
entirely  clear,  but  most  of  the  error  occurs  during 


Table  3.  Errors  in  frost  depth  calculations  using  linear  shape  functions.* 


Analytical  solution 
t lapsed  (modified  Berqgren 

time  equation)  Difference  between  analytical  solution  and  FEM  solution  (analytical  -  FEM). 

(h) _ (err) _ (err) _ 


Node  spacing  (cm) 

0.5 

2.0 

2.0 

1.0 

2.0 

5.0 

10.0 

Time  increment  (h) 

0.02 

0.1 

0.25 

0.5 

0.5 

0.5 

0.5 

1 

2.28 

-0.06 

-0.45 

+0.14 

1.00 

O.tS 

-0.02 

-1.57 

2 

3.23 

-0.06 

-0.77 

-0.01 

0.94 

0.73 

-1.37 

-4.46 

3 

3.96 

-0.06 

-0.67 

-0.13 

0.94 

0.56 

-1.32 

-6.04 

4 

4.57 

-0.06 

-0.75 

-0.14 

0.96 

0.45 

-1.27 

-5.78 

5 

5.11 

-0.07 

-0.89 

-0.21 

0.93 

0.45 

-1.29 

-5.59 

12 

7.91 

-0.10 

-0.82 

-0.30 

0.91 

0.23 

-2.10 

-5.16 

24 

11.19 

-0.10 

-0.89 

-0.41 

0.96 

0.23 

-2.08 

-5.87 

48 

15.82 

-0.12 

-0.90 

-0.42 

0.92 

0.10 

-2.23 

-5.84 

72 

19.38 

-0.13 

-0.92 

-0.50 

0.94 

0.11 

-2.24 

-5.89 

♦Soil  properties: 

Volumetric  heat  capacity  =  1.007  cal/cmJ  °C 
Thermal  conductivity  =  16.05  cal/cm  h  °C 
Latent  heat  of  fusion  =  28.96  cal/cm1 


Initial  conditions: 

Semi-infinite  body  isothermal  at  0°C 
Boundary  condition: 

Temperature  -  -5°C  at  surface  for  time  >  0 


freezing  of  the  uppermost  element.  One  factor  con¬ 
tributing  to  the  inaccuracy  of  the  FEM  solution  is 
that  the  forced  upper  boundary  conditions  do  not 
precisely  represent  the  analytical  problem.  The  an¬ 
alytical  problem  is  one  where  the  semi-infinite  body 
is  initially  at  uniform  temperature  and  for  succeeding 
times  a  different  temperature  is  applied  at  the  upper 
surface.  In  contrast,  the  FEM  solutions  can  only 
approximate  this  step  change,  because  a  smooth 
gradient,  defined  by  the  shape  function,  exists  between 
the  temperature  at  the  surface  node  and  that  at  the 
other  boundary  node  of  the  uppermost  element.  The 
freezing  algorithm  requires  repeated  application  of  this 
approximation  if  the  element  takes  more  than  one  time 
increment  to  freeze. 

Despite  some  error  near  the  top  of  the  profile,  dif¬ 
ferences  of  only  a  few  percent  generally  occur  after 
several  hours  of  simulation.  Similar  studies  were 
conducted  using  second  order  shape  functions,  and 
behavior  with  respect  to  various  combinations  of  the 
time  and  depth  increments  was  similar;  the  accuracies 
of  the  solutions  were  not  significantly  improved  over 
those  obtained  with  linear  shape  functions. 

Moisture  Flux 

There  are  no  known  analytical  solutions  for  the 
Richard’s  equation  for  unfrozen  or  freezing  soils. 
Therefore,  to  verify  the  finite  element  representation 
of  the  equation,  above-freezing  isothermal  conditions 
were  imposed  and  comparisons  were  made  with  an 
approximate  solution  suggested  by  Philip  (1957)  and 
detailed  bv  Kirkham  and  Powers  (1972).  Guymon  and 


Luthin  (1974)  state  that  rough  computations  by  Philip’s 
method  using  the  same  soil-water  parameters  indicate 
that  the  wetting  curves  in  their  report  were  “approx¬ 
imately  correct."  Two  of  the  curves  from  Guymon 
and  Luthin  (1974)  are  shown  in  Figure  6  as  are  two 
curves  from  the  present  computer  model.  Agreement 
between  the  two  models  is  very  good,  and  data  from 
these  and  other  simulations  indicate  that  the  computer 
model  for  moisture  flux  under  isothermal  conditions 
approximates  the  test  solution  reasonably  well.  It  was 
also  concluded  that  the  finite  element  solution  was  not 
particularly  sensitive  to  parameter  update  frequency 
or  element  length. 

Several  test  solutions  were  conducted  to  determine 
the  sensitivity  of  the  finite  element  model  to  parameter 
errors.  Parameters  evaluated  included  hydraulic  con¬ 
ductivity,  porosity  and  Gardner’s  constants  Ak  and  /4W. 
Results,  which  were  developed  from  an  early  version  of 
the  model  using  linear  shape  factors,  are  shown  in  Figure 
7.  The  ordinate  in  Figure  7  is  the  average  difference  in 
computed  pressures  for  the  "base”  case  and  the  altered 
case.  The  average  pressure  is  determined  from  the 
pressures  at  depths  of  1 0,  20,  30  and  40  cm.  The  data 
indicate  that  solution  accuracy  is  highly  dependent  on 
the  hydraulic  conductivity  and  porosity.  Data  also  in¬ 
dicate  that  a  change  in  the  hydraulic  conductivity  has 
the  opposite  effect  of  a  change  in  the  porosity  and  that 
changes  in  A  w  and  A  k  also  have  opposing  effects.  Ad¬ 
ditional  discussion  on  the  effect  of  Gardner’s  coefficients 
and  the  hydraulic  conductivity  is  in  Frost  Heave  of 
Homogeneous  Laboratory  Samples, 
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Water  Content  tern s/em*) 


Figure  6.  Comparison  of  FEM  solutions  of  Richard's  equation. 
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Figure  7.  Sensitivity  of  parameter  errors. 


Numerical  Dispersion 

The  heat  flow  and  moisture  flow  equations  are  based 
on  the  assumption  that  the  primary  flux  is  by  conduction 
of  heat  and  moisture.  Significant  numerical  dispersion 
may  be  introduced  into  the  solutions  of  these  nonsym- 
metrical  equations  if  convective  effects  dominate  or  are 
overly  significant.  In  the  case  of  Richard’s  equation, 
these  effects  are  determined  by  the  dimensionless  ratio, 
K*  Ax/K,  where  K  is  the  unsaturated  hydraulic  con¬ 
ductivity,  K*  is  defined  by  eq  6a,  and  Ax  is  the  element 
length.  If  this  ratio  reaches  a  value  near  unity,  numerical 
dispersion  may  become  significant.  For  most  actual 
conditions,  it  is  not  anticipated  that  numerical  dispersion 
caused  by  convective  effects  will  be  a  problem.  However, 
in  our  sensitivity  studies  this  ratio  did  at  times  reach  a 


critical  state  at  one  or  more  nodal  points.  The  ratio  was 
calculated  for  each  nodal  point  at  each  time  step  and  if 
the  ratio  reached  or  exceeded  unity  the  computer  solution 
was  terminated. 

For  the  heat  flux  equation,  the  appropriate  critical 
ratio  is  Cw  V  Ax/Ky  where  Cw  is  the  volumetric  heat 
capacity  of  water,  V  the  velocity  flux,  and  Ky  the  thermal 
conductivity.  As  in  the  moisture  flux  equation,  numerical 
dispersion  may  become  significant  if  the  ratio  approaches 
unity.  This  ratio  is  also  computed  at  each  nodal  point 
at  each  time  increment.  If  the  ratio  reaches  or  exceeds 
unity,  the  solution  is  terminated. 

For  a  typical  soil,  Cw  equals  about  1  cal/cm3  °C  and 
Ky  =  5.0x  1  O'3  cal/cm  s  °C.  An  element  length  of  2  cm 
is  typical  in  model  trials.  If  unity  is  taken  as  the  critical 
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Table  4.  Comparison  of  calculated  and  experimental  frost  penetration  and  frost  heave. 
All  frost  depths  are  from  the  original  surface.  Soil  Is  New  Hampshire  silt  with  an  initial  density  of 
1 .73  g  cm'3  and  an  initial  moisture  content  of  3 1 .3%  by  volume. 

Calculated  results 


At  =  1.0 


■op  3 

A  k  =  7.0X10'°  cm'3 


At  =  1.0  A 

A  k  =  7.0X1 0'7  cm'3 
<31  3=  1.83x10"*  cm/h  K3]  3  =  1.7X10'3  cm/h 


At  =  0.75/1 
Ak  =  7.0X1 0'7  cm'3 
*  ,3 
K31  3  =  1.75X10  3  cm/h 


Experimental  Results 


Convection  In 
both  equations 
(Simulation  1) 


Convection  In 
both  equations 
(Simulation  2) 


Convection  In 
both  equations 
(Simulation  3) 


No  heat  transfer 
by  convection 
(Simulation  4) 


Neither  heat 
nor  mass 
transfer  by 
convection 
(Simulation  S) 


Elapsed 

Frost 

Frost 

Frost 

Frost 

Frost 

Frost 

Frost 

Frost 

Frost 

Frost 

Frost 

Frost 

time 

depth 

heave 

depth 

heave 

depth 

heave 

depth 

heave 

depth 

heave 

depth 

heave 

(days) 

(cm) 

(cm) 

(cm) 

(cm) 

S1H33J1 

want m 

(cm) 

(cm) 

(cm) 

(cm) 

(cm) 

(cm) 

1.0 

2.1 

0.6 

6.0 

0.3 

6.0 

0.5 

6.0 

0.6 

6.0 

.6 

6.0 

.6 

2.0 

- 

- 

8.2 

0.5 

8.0 

0.9 

8.6 

1.0 

8.6 

1.0 

8.2 

1.1 

3.0 

6.1 

1.3 

10.8 

0.6 

10.0 

1.3 

10.0 

1.6 

10.0 

1.6 

10.0 

1.6 

4.0 

7.1 

1.5 

11.6 

0.9 

10.0 

2.1 

12.0 

2.0 

12.0 

2.0 

10.0 

IS 

5.0 

7.7 

1.7 

13.7 

0.9 

10.0 

3.0 

12.0 

2.8 

12.0 

STB 

10.0 

3.6 

ratio,  then  a  velocity  of  2.5x  10  3  cm/s  is  critical.  Silts 
and  slightly  coarser-grained  soils  may  have  saturated 
hydraulic  conductivity  (permeability)  values  in  this 
range.  This  ratio  is  also  computed  for  each  nodal  point 
at  each  time  step.  If  the  ratio  reaches  or  exceeds  unity, 
the  computer  solution  is  terminated.  Guymon  and 
Hrodmadka  (in  prep.)  are  studying  a  slightly  different 
formulation  which  virtually  eliminates  this  problem. 

Frost  Heave  of  Homogeneous  Laboratory  Samples 

Initial  solutions  to  problems  with  arbitrarily  chosen 
initial  and  boundary  conditions  indicated  that  the  frost 
heave  algorithm  behaved  well  qualitatively.  Therefore, 
comparison  was  made  between  frost  depths  and  frost 
heave  computed  by  the  model  and  values  measured  in 
controlled  laboratory  experiments.  The  measured 
data  were  from  a  standard  CRREL  frost-susceptibility 
test  on  Manchester  silt  (generally  called  New  Hamp¬ 
shire  silt  in  the  literature).  Initial  conditions  for  the 
model  were  measured  conditions  at  the  beginning  of 
the  test,  and  the  time  dependent  temperature  bound¬ 
ary  conditions  for  the  model  were  those  measured  dur¬ 
ing  the  test. 

Since  unsaturated  hydraulic  conductivities  at  different 
moisture  contents  had  not  been  measured  for  New 
Hampshire  silt,  the  value  of  Ak  (used  in  Gardner's 
relationship)  was  estimated.  Accuracy  of  the  estimate 
was  determined  from  a  comparison  of  calculated  and 
experimental  frost  heave  and  frost  penetration.  Table 


4  contains  a  comparison  of  experimental  and  calculated 
results.  The  same  spatial  discretization,  2-cm  element 
length,  was  used  for  all  simulations  in  Table  4. 

Comparison  of  simulations  1  and  2  shows  the  effect 
of  changing  Gardner’s  soil  permeability  coefficient  A  k . 
The  saturated  hydraulic  conductivity  (permeability)  of 
this  soil  was  4.0x  10'J  cm/h  at  a  porosity  of  36.2%. 

Using  these  data  in  eq  4,  hydraulic  conductivities  of 
1 .83x  10'4  cm/h  at  a  moisture  content  of  31 .3%  by 
volume  (300  cm  of  water  tension)  and  of  6.65x  10"6 
cm/h  at  a  moisture  content  of  6.8%  by  volume  (950 
cm  of  water  tension)  are  computed.  The  value  of  Ak 
used  in  simulations  2-5  corresponded  to  a  hydraulic 
conductivity  of  1.75x1 0"3  cm/h  at  the  31.3%  moisture 
content  and  to  6.64x  10'5  cm/h  at  the  6.8%  moisture 
content.  Therefore,  between  simulation  1  and  2  un¬ 
saturated  hydraulic  conductivities  differ  by  roughly  an 
order  of  magnitude.  In  simulation  1 ,  the  lower  hydraul¬ 
ic  conductivities  resulted  in  the  complete  cessation  of 
heave  after  four  days  (while  the  freezing  front  advanced 
steadily);  in  simulation  2,  a  massive  ice  lens  grew  after 
three  days  and  the  freezing  front  ceased  to  advance  be¬ 
cause  the  rate  of  heat  extraction  was  just  adequate  to 
freeze  the  inflowing  moisture. 

Comparison  of  simulations  2  and  3  illustrates  the 
effect  of  time  discretization.  In  going  from  one  time 
increment  size  to  another,  the  rate  of  freezing  is  only 
slightly  altered  in  the  second  day  of  simulation.  However, 
when  the  rate  of  heat  extraction  is  slow  and  water  influx 


Figure  8.  Test  cylinder  to  evaluate  effects  of  surcharge  on  frost  heaving. 
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is  significant,  the  freezing  rate  is  strongly  coupled  with 
the  existing  spatial  discretization,  so  that  relatively  large 
differences  exist  between  the  two  simulations  for  days 
4  and  5. 

Simulations  3,  4  and  5  illustrate  the  magnitude  of 
the  contribution  of  the  convective  terms.  By  eliminating 
the  convective  contribution  to  heat  transport  (simulation 
4)  we  see,  at  most,  only  a  small  change  in  calculated 
temperatures,  frost  depths,  and  the  amount  of  heave. 
Elimination  of  the  convection  term  in  the  moisture 
flow  equation  (simulation  S)  has  a  slightly  larger  effect 
on  the  first  few  days  of  simulated  results  and  a  dramatic 
effect  of  misleading  magnitude  on  the  results  beginning 
with  day  4,  which  were  dominated  by  the  influence  of 
spatial  discretization.  The  early  days  of  simulation  show 
slower  frost  penetration  and  greater  frost  heaving  in 
simulation  5  than  in  simulation  4  because  the  convective 
mass  transport  term  includes  the  effect  of  gravity;  its 
elimination  allows  more  water  to  be  drawn  to  the  freezing 
front. 

A  test  cylinder  (Fig.  8)  was  constructed  to  develop 
data  to  provide  a  preliminary  verification  of  the  overburden 
and/or  surcharge  portion  of  the  frost  heave  algorithm. 

The  test  cylinder  was  about  1  m  longx  1  S-cm  inside 
diameter.  The  upper  1 5  cm  of  the  cylinder  was  placed 
inside  a  freezing  cabinet.  It  is  removable  from  the 
lower  portion.  Loose  insulation  was  placed  around  the 
upper  15  cm  segment  and  only  the  upper  surface  was 
exposed  to  freezing  temperatures,  thus  allowing  one¬ 
dimensional  heat  flux  from  the  upper  portion  of  the 
cylinder.  Tensiometers  and  thermocouples  were  in¬ 
stalled  within  the  soil  to  measure  soil  water  tension  and 
temperatures,  respectively.  A  “bubble  tube”  reservoir 
was  used  to  control  the  free-water  level  in  the  sample. 

Berg  et  al.  (in  prep.)  discuss  the  test  cylinder  and 
results  from  three  tests  which  were  conducted  with 
the  cylinder.  Figures  9  and  10  show  results  from  two 
tests  in  which  a  surcharge  equivalent  to  35  cm  of  water 
was  used.  Fairbanks  silt  compacted  to  a  density  of 
approximately  1 .55  g/cm3  was  used  in  all  tests.  Samples 
were  placed  and  compacted  slightly  below  the  optimum 
moisture  content.  Prior  to  initiating  freezing,  the  water 
table  was  raised  to  the  surface  of  the  sample  until  free 
water  appeared.  It  was  then  lowered  to  a  depth  of  45 
cm  for  several  days  prior  to  freezing.  It  was  held  at  the 
45-cm  depth  during  the  test. 

Since  the  initial  temperature  and  moisture  conditions 
were  similar  in  both  tests,  a  series  of  computer  simulations 
were  made  and  compared  with  measured  data  after  10 
days  of  freezing;  results  are  shown  in  Table  5.  Various 
combinations  of  three  parameters  affecting  the  quantity 
of  moisture  flux  and  thus  the  magnitude  of  frost  heaving 
were  used  in  the  simulations.  The  three  parameters 
varied  were  Gardner’s  hydraulic  conductivity  coefficient 


4k  Gardner’s  moisture  coefficient  A  w  and  the  per¬ 
meability  or  saturated  hydraulic  conductivity  (eq  4-7). 

Simulation  1  used  values  for  the  three  parameters 
which  were  measured  on  samples  of  Fairbanks  silt  at 
approximately  the  same  density  as  the  samples  in  the 
freezing  tests.  With  these  values  neither  the  frost  depth 
nor  the  amount  of  frost  heave  was  accurately  computed. 
Several  additional  simulations  were  made  by  changing 
the  value  of  one  or  more  parameters. 

In  simulations  2-7,  only  one  parameter  yi r  simulation 
differed  from  values  used  in  simulation  1.  In  simulations 
8-16,  more  than  one  of  the  parameters  were  varied  in 
each  trial.  Results  from  simulations  4  and  7  indicate 
that  increasing  the  permeability  enhanced  both  the 
frost  depth  and  the  frost  heave  after  10  days.  The 
amount  of  frost  heave,  however,  was  only  about  one- 
third  of  that  observed  in  the  laboratory  tests.  Results 
from  simulation  15  agreed  most  closely  with  measured 
values  of  frost  heave,  but  computed  and  measured 
frost  penetration  depths  do  not  agree  as  well  as  in  some 
other  trials. 

Figure  1 1  illustrates  computed  and  measured  values 
of  frost  heave  and  frost  penetration  as  functions  of  time. 
Frost  heave  is  computed  reasonably  well  throughout  the 
entire  duration  of  the  test.  The  upper  set  of  ’’calculated" 
frost  depth  points  represents  the  lower  limit  of  the 
solidly  frozen  zone;  i.e.  computed  temperatures  above 
this  level  are  below  the  freezing  point  depression  tem¬ 
perature  (-0.1  °C)  in  the  simulations  for  this  soil  or, 
stated  in  another  way,  sufficient  excess  degrees  (eq  27) 
have  been  accumulated  to  freeze  the  soil  to  this  depth. 
The  lower  set  of  "calculated”  frost  depth  points  forms 
the  lower  limit  of  the  0°C  isotherm.  In  other  words, 
the  area  between  the  two  sets  of  “computed”  frost 
depth  points  is  a  freezing  zone  as  determined  by  the 
model.  Calculated  temperatures  in  the  freezing  zone 
are  between  0°C  and  -0.1  °C. 

Electrical  resistivity  gages  (Atkins  1976)  were  in¬ 
stalled  in  the  sample  used  in  test  2.  These  devices 
operate  on  the  principle  that  the  electrical  resistivity 
of  ice  is  much  greater  than  that  of  water.  In  test  2, 

copper  wires  were  flattened  and  installed  1  cm  below  the 
sample  surface  and  at  2-cm  intervals  to  a  depth  of  1 5  cm. 

By  measuring  the  electrical  resistance  between  each 
pair  of  r-ectrodes  and  plotting  the  data  as  a  function  of 
time,  it  was  possible  to  determine  the  time  when  about 
one-half  of  the  soil  moisture  was  frozen  and  the  time 
when  most  of  the  remaining  moisture  froze.  (This  soil 
has  a  significant  amount  of  unfrozen  water  at  sub¬ 
freezing  temperatures  near  0°C,  according  to  Anderson 
and  Tice  1972.)  Figure  12  contains  plots  of  the  0°C 
isotherm  vs  time  as  well  as  lines  showing  when  approxi¬ 
mately  50%  and  100%  of  the  soil  water  was  frozen  at 
various  depths. 
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Figure  9.  Results  from  test  1  on  Fairbanks  silt  with  surcharge  equivalent  to  35  cm  of  water. 


MjfeQ  0 


a.  Temperatures,  water  absorption,  frost  heave  and  frost  b.  Soil  moisture  stress  at  various  depths, 

depth  plots. 


Fro*t  Dtp 


r 


Table  S.  Frost  penetration  and  frost  heave  after  10  days  of 
laboratory  freezing. 


Simu¬ 

lation 

»k* 

t 

Permeability 

(cm/h) 

Frost 

depth  (cm) 

Frost 

heave  (cm) 

1 

5E-5** 

3E-9 

5E-S 

9.69 

0.06 

2 

5E-4 

3E-9 

5E-5 

9.39 

0.00 

3 

5E-5 

3E-8 

5E-5 

10.00 

0.15 

4 

5E-5 

3E-9 

5E-4 

10.00 

0.33 

5 

5E-5 

3E-7 

5E-5 

11.20 

0.09 

6 

5E-6 

3E-9 

5E-S 

9.98 

0.17 

7 

5E-5 

3E-9 

5E-3 

10.50 

0.66 

8 

5E-4 

3E-9 

5E-3 

10.00 

0.3S 

9 

5E-3 

3E-9 

5E-3 

9  69 

0.06 

10 

SE-4 

3E-10 

5E-3 

7.78 

0.00 

12 

5E-6 

3E-9 

5E-3 

9.99 

1.24 

13 

5E-5 

3E-8 

5E-3 

10.30 

0.98 

14 

5E-6 

3E-8 

5E-3 

10.50 

0.89 

15 

5E-6 

3E-9 

5E-2 

8.5 

2.28 

16 

5E-6 

3E-9 

1.1E-1 

7.9 

3.69 

- 

Laboratory  test  1 

12.0 

2.75 

- 

Laboratory  test  2 

12.5 

2. OS 

♦Gardner’s  hydraulic  conductivity  coefficient. 
fGardner’s  moisture  coefficient. 

»»5E-5  =  5X10'S  etc. 


Figure  1 2.  Temperature  and  electrical  resistivity  gage  data  for  test  2. 


Figure  1 1 .  Computed  and  measured  frost 
heave  and  frost  penetration  on  Fairbanks  silt 
with  surcharge  of  35  cm  of  water. 
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a.  Temperatures,  water  absorption,  frost  heave  and  frost  depth  plots. 

Figure  13.  Results  from  test  3  on  Fairbanks  silt  with  surcharge  equiva¬ 
lent  to  350  cm  of  water. 


The  freezing  zone,  i.e.  the  distance  between  the 
0°C  isotherm  and  the  solidly  frozen  line,  is  not  nearly 
as  wide  as  the  computed  freezing  zone  in  Figure  1 1 . 

The  width  of  the  zone  is  not  constant,  but  increases 
with  time.  Until  the  solidly  frozen  soil  penetrated  to 
a  depth  at  about  7  cm,  it  corresponded  to  a  temperature 
of  about  -0.3°C;  thereafter  the  boundary  of  the  solidly 
frozen  soil  corresponded  with  progressively  lower  tem¬ 
peratures.  On  3  May  the  solidly  frozen  boundary  was 
at  13  cm  which  conformed  to  a  temperature  of  -2.2°C. 
An  unintentional,  partial  freeze-thaw  cycle  occurred 
on  26  April  (Fig.  11)  when  die  sample  thawed  to  a  depth 
of  about  5  cm  and  the  remainder  of  it  warmed  sub¬ 
stantially.  Whether  this  cycle  affected  the  amount  of 
unfrozen  water  when  it  was  refrozen  is  unknown  at 
this  time. 

Data  from  tests  1  and  2  and  from  the  computer 


simulations  in  Table  5  indicate  that  the  model  can  ac¬ 
curately  predict  the  amount  of  frost  heave  and  frost 
penetration  in  short-duration  laboratory  samples.  Data 
from  Table  5  again  illustrate,  however,  that  heave  and 
frost  penetration  are  quite  sensitive  to  the  hydraulic 
properties  used  in  the  solutions. 

A  third  test  was  conducted  in  this  series  to  determine 
the  validity  of  eq  32.  A  surcharge  equivalent  to  a 
pressure  of  350  cm  of  water  was  applied  to  the  soil 
column  prior  to  initiation  of  freezing.  If  eq  32 
adequately  represented  effects  of  overburden  and 
surcharge,  the  soil  water  tension  measurements  would 
have  been  lower,  i.e.  less  negative,  during  this  test  than 
during  the  previous  two  tests. 

Figure  1 3,  which  contains  data  observed  during  the 
third  test,  shows  that  the  measured  soil  water  tension, 
values  are  not  significantly  different  from  those  meas- 
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b.  Soil  moisture  stress  at  various  depths. 

Figure  13.  (cont’d) 

ured  in  the  other  two  tests  with  a  surcharge  pressure  CONCLUSIONS 

of  35  cm  of  water  (Fig.  9  and  10).  Data  from  Figure 

10  indicated  that  the  maximum  soil  water  tension  ob-  A  one-dimensional  mathematical  model  coupling 

served  as  the  freezing  front  approached  a  tensiometer  at  simultaneous  heat  and  moisture  fluxes  in  freezing  soil- 

a  depth  of  15  cm  was  about -650  cm  of  water.  Using  water  systems  has  been  developed.  The  finite  element 

this  value  of  Pw  in  eq  32  we  calculate  Pc  =  -1622  cm  method  was  used  to  solve  the  problem,  and  others  (e.g. 

of  water  for  test  2.  If  we  assume  PQ  does  not  change  Taylor  and  Luthin  1 976  and  1 978,  Harlan  1 972  and 

due  to  the  added  surcharge,  we  anticipate  that  a  maxi-  1973,  Dempsey  1978  and  Outcalt  1976)  have  solved 

mum  tensiometer  reading  of  about  -490  cm  of  water  the  one-dimensional  coupled  heat  and  moisture  flow 

will  occur  at  the  freezing  front  in  test  3.  Unfortunately,  equations  using  finite  difference  procedures, 
a  tensiometer  was  not  located  in  the  freezing  zone  in  This  wodel  has  successfully  been  used  to  simulate 

test  3,  but  tensions  at  the  1 8-cm  depth  and  below  did  frost  penetration  and  frost  heave  which  were  observed 

not  show  a  significant  change  from  test  2.  These  results  in  relatively  homogeneous  laboratory  test  specimens, 

do  not  completely  eliminate  the  possibility  that  eq  32  The  amount  of  water  flux,  and  resultant  frost  heave, 

is  valid,  however,  because  soil  moisture  tension  data  computed  from  the  model,  is  very  sensitive  to  the 

were  not  measured  in  the  immediate  vicinity  of  the  hydraulic  properties  of  the  soil.  This  agrees  qualitatively 

freezing  zone  in  all  three  tests.  It  is  possible  that  sig-  with  field  observations;  i.e.  certain  soils,  such  as  silts, 

nificant  changes  in  tension  values  in  this  area  were  not  have  optimum  hydraulic  properties  to  cause  severe 

reflected  in  pressure  gradients  a  few  centimeters  away.  frost  heave. 

Berg  et  al.  (in  prep.)  present  a  thorough  discussion  A  work  plan  was  developed  for  correlating  results  of 

of  results  from  all  three  tests  which  were  conducted  in  various  laboratory  test  methods  for  assessing  frost  heave 

this  series.  susceptibility  of  soils  and  the  amount  of  heave  observed 
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at  field  study  sites.  This  mathematical  model  can  be 
used  to  help  correlate  the  field  and  laboratory  data. 
Appendix  A  contains  the  work  plan. 

Field  instrumentation  necessary  for  correlating 
laboratory  and  field  observations  of  frost  heave  would 
be  nearly  the  same  as  that  required  for  investigations 
of  the  thaw-weakening  of  subgrade  soils  and  granular 
unbound  base  courses.  A  proposed  investigational 
procedure  for  laboratory  and  field  studies  related  to 
thaw  weakening  is  presented  in  Appendix  B. 

Possibly  the  most  important,  but  as  yet  unsubstan¬ 
tiated,  application  of  the  mathematical  model  is  to 
predict  the  amount  of  frost  heave  which  may  occur  in 
a  given  pavement  system  at  a  given  geographical  location. 
Since  it  is  based  entirely  upon  principles  of  heat  and 
moisture  flux,  the  model  does  not  require,  or  use,  data 
from  any  of  the  laboratory  frost-susceptibility  test 
methods.  If  the  mathematical  model  can  be  used  to 
estimate  the  likely  range  of  frost  heave  of  candidate 
pavement  profiles,  we  have  made  a  significant  improve¬ 
ment  in  the  procedure  for  designing  roadway  and  air¬ 
field  pavements  in  seasonal  frost  areas. 


RECOMMENDED  STUDIES  TO  REFINE  THE  MODEL 

In  September  1978,  the  Corps  of  Engineers,  the 
Federal  Highway  Administration  and  the  Federal  Avi¬ 
ation  Administration  commenced  a  study  to  correlate 
laboratory  and  field  measurements  of  frost-susceptibility 
and  frost  heave  and  laboratory  and  field  studies  of  the 
extent  and  duration  of  thaw  weakening.  The  model 
developed  in  the  present  study  will  be  refined  and 
applied  in  the  recently  initiated  study. 

Necessary  refinements  to  the  mathematical  model 
can  be  divided  into  two  broad  categories:  1 )  to  improve 
computational  efficiency  and  2)  to  improve  the  accuracy 
of  estimates  of  frost  penetration  and  frost  heave.  Im¬ 
provements  in  the  computational  efficiency  should 
include  more  efficient  storage  of  data  within  the 
matrices,  less  time-consuming  manipulation  of  the 
matrices  and  improved  formulation  of  certain  aspects 
of  the  problem  to  allow  using  smaller  computers  or 
decreasing  solution  times. 

Improving  the  accuracy  of  estimating  frost  penetration 
and  frost  heave  requires  a  more  complete  understanding 
and/or  representation  of  thermal,  hydrological  and 
mechanical  processes  within  the  freezing  zone.  Specif¬ 
ically,  algorithms  which  incorporate  the  effects  of  sur¬ 
charge  and  overburden  are  necessary.  Expressions  for 
the  relationships  between  soil  moisture  tension  and 
moisture  content  and  soil  moisture  tension  and  un¬ 
saturated  hydraulic  conductivity  are  used  in  the  com¬ 
putations.  We  must  develop  expressions  which  suitably 


relate  these  parameters,  and  it  may  also  be  necessary 
to  incorporate  effects  of  hysteresis  on  these  parameters. 

In  the  recently  initiated  joint  study  it  is  planned  to 
construct  a  dual  gamma  beam  device  which  can  non¬ 
destructive^  scan  through  a  soil  column  similar  to  the 
one  shown  in  Figure  8.  This  device  will  permit  simul¬ 
taneous  measurement  of  changes  in  moisture  content 
and  density  as  the  soil  column  freezes.  In  addition, 
relationships  between  soil  moisture  tension  and  moisture 
content  and  between  soil  moisture  tension  and  un¬ 
saturated  hydraulic  conductivity  can  be  determined 
directly  from  data  obtained  with  the  soil  column  in¬ 
strumentation  and  the  dual  gamma  system. 

Laboratory  tests  involving  freeze-thaw  cycles  and 
layered  soil  systems  may  be  required  to  more  thoroughly 
validate  the  model  prior  to  predicting  frost  penetration 
and  frost  heave  from  actual  road  and  airfield  pavements. 
These  more  complex  problems  have  not  been  addressed 
in  the  current  study,  and  although  no  major  problems 
are  anticipated  when  extending  the  model  to  these 
conditions,  minor  complications  will  inevitably  arise. 

The  final  test  of  the  model  will  be  to  predict  frost 
heave  and  frost  penetration  of  instrumented  field  test 
sections.  If  the  model  passes  these  tests,  it  can  confidently 
be  applied  to  assessing  frost  heave  in  various  pavement 
profiles  under  various  environmental  conditions. 

This  one-dimensional  model  cannot  be  used  to  pre¬ 
dict  differential  frost  heave.  Perhaps  the  model  should 
be  used  in  conjunction  with  another  model,  which  may 
be  statistically  based,  which  describes  spatial  variations 
of  the  important  parameters.  This  problem  will  be 
studied  during  the  recently  initiated  studies. 
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APPENDIX  A:  WORK  PLAN,  STAFFING  AND 
INSTRUMENTATION  REQUIREMENTS  FOR  COR¬ 
RELATING  RESULTS  OF  LABORATORY  FROST 
SUSCEPTIBILITY  TESTS  WITH  FIELD  PERFOR¬ 
MANCE 

Excerpted  from  July  1976  Interim  Report  by  R.L. 

Berg  and  T.C.  Johnson. 

INTRODUCTION 

In  July  1975,  the  Federal  Highway  Administration 
(FHWA),  the  Federal  Aviation  Administration  (FAA) 
and  CRREL  started  a  jointly  funded  project  to  develop 
a  mathematical  model  of  the  frost  heaving  process.  The 
work  is  being  accomplished  by  CRREL. 

Prior  to  conducting  field  studies,  it  was  necessary  to 
develop  a  numerical  model  for  predicting  frost  heave 
in  roadway  and  airfield  pavements.  Otherwise  there 
would  have  been  no  effective  way  to  normalize  other 
critical  parameters  at  the  various  test  sites  to  permit 
comparative  evaluation  of  the  different  laboratory 
methods  of  determining  soil  frost-susceptibility.  The 
model  will  also  permit  important  parameters  to  be 
varied  separately  or  simultaneously  to  assess  their  in¬ 
fluence  on  frost  heaving  and  should  be  a  valuable  ana¬ 
lytical  tool  for  use  by  designers  in  predicting  frost 
heave. 

Development  of  the  mathematical  model  is  the  ini¬ 
tial  step  in  a  broader  study  which  will  ultimately  in¬ 
volve  State  Highway  Departments,  State  Departments 
of  Transportation,  and  airfield  construction  agencies 
in  a  cooperative  undertaking  to  measure  frost  penetra¬ 
tion,  frost  heave  and  other  important  variables  in  full- 
scale  test  sections.  When  the  field  tests  are  conducted, 
various  laboratory  methods  for  assessing  the  frost  sus¬ 
ceptibility  of  the  soils  will  be  applied  to  determine 
which  method  gives  the  most  accurate  or  consistent 
index  for  predicting  the  magnitude  of  frost  heave. 

The  remainder  of  this  appendix  discusses  a  suggested 
work  plan  for  the  field  studies,  including  the  composi¬ 
tion  of  a  technical  team  to  conduct  tl  studies,  general 
criteria  for  site  selection,  and  suggested  instrumentation 
for  each  test  section. 


FIELD  TEST  SITES 

The  validity  of  this  frost  model  will  be  established 
by  full-scale  test  sections  installed  at  selected  locations 
in  several  of  the  northern  states.  It  would  be  desirable 
to  provide  centralized  direction  for  the  program  by 


making  a  single  organization  responsible  for  coordinating 
the  overall  study,  providing  assistance  on  instrumentation 
and  observations,  directing  and  coordinating  the  related 
laboratory  studies,  and  interpreting  and  analyzing 
results.  The  organization  should  have  a  team  of  spe¬ 
cialists  with  expertise  in  the  following  areas:  1 )  labora¬ 
tory  soil  tests  for  measurement  of  frost-susceptibility 
and  of  related  soil  properties,  2)  instrumentation  for 
laboratory  and  field  measurement  of  important  thermal 
and  hydrological  parameters,  3)  collection,  storage  and 
management  of  data  from  field  tests,  4)  theories  of  the 
frost  heaving  process,  5)  physics  and  chemistry  of  soil- 
water  processes,  6)  meteorology  and  computational 
methods  for  assessing  the  surface  energy  balance  on 
pavement  surface,  7)  soil  mechanics,  8)  surface  and 
subsurface  drainage  of  pavement  systems,  9)  pavement 
design  and  construction  in  seasonal  frost  areas,  and 
10)  computer  programing.  In  addition,  the  organiza¬ 
tion  coordinating  the  study  should  have  the  facilities 
and  resources  necessary  to  process  the  large  volumes  of 
data  that  will  be  generated  from  the  laboratory  and 
field  studies. 

Selection  of  appropriate  existing  pavements  for  the 
field  test  sections  and  a  suitable  choice  of  instrumenta¬ 
tion  are  critical  tasks  in  the  proposed  research  to  assure 
the  highest  probability  of  correlation  of  data  from  the 
various  field  sites,  the  laboratory  test  data,  and  the 
mathematical  model.  To  minimize  the  amount  of  in¬ 
strumentation  required,  and  to  keep  the  volume  of 
data  generated  within  manageable  limits,  it  is  extremely 
important  that  sites  be  chosen  where  there  is  essenti¬ 
ally  no  lateral  or  longitudinal  flow  of  heat  or  moisture. 
Sites  of  this  type  would  lend  themselves  to  a  better 
simulation  by  both  the  laboratory  tests  and  the  math¬ 
ematical  model  and  would  require  less  costly  instrumen¬ 
tation  and  data  collection.  Table  A1  lists  factors  that 
should  be  considered  in  the  site  selection  process. 

Figures  A1  and  A2  illustrate  the  observation  points 
and  the  instrumentation  to  be  installed  at  a  typical  test 
site  in  a  highway  pavement.  In  an  airfield  pavement 
the  instrumentation  would  be  similar  but  the  test  pit 
could  be  located  1  5-20  ft  (4.6-6. 1  m)  from  a  runway 
edge.  The  amount  and  location  of  instrumentation 
are  idealized  in  the  two  figures,  and  in  a  definitive  study 
should  be  tailored  to  specific  site  conditions,  availability 
of  equipment  and  personnel,  and  funding.  These  same 
factors  may  also  strongly  influence  the  type  and  frequency 
of  observations  to  be  made.  Details  for  instrumentation 
at  a  specific  site  should  be  developed  through  discussions 
between  the  sponsors,  the  involved  highway  or  airport 
operators,  and  the  coordinating  organization.  The  types 
of  instrumentation  and  observations  outlined  below  are 
suggested  as  an  optimum  system  that  will  provide 
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Table  A1.  Considerations  in  site  selection. 


1 .  Little  or  no  transverse  slope  except  crown  in  pavement. 

2.  Little  or  no  longitudinal  slope. 

3.  Pavement  section  not  excessively  layered,  i.e.  few  materials  of  widely  differing 
properties. 

4.  Pavement  surface  not  severely  cracked. 

5.  Stable  water  table. 

6.  Occurrence  of  frost  heave  of  the  pavement  surface. 

7.  Pavement  kept  cleared  of  snow  in  winter. 

8.  Uniform  subsurface  soil  and  moisture  conditions. 

9.  Uniform  exposure  to  solar  radiation  and  wind. 

10.  Minimum  of  100-ft-long  (30.5  m)  pavement  segment. 

11.  Within  a  reasonable  distance  of  office  responsible  for  observations  and  servicing 
of  instrumentation. 

12.  Other  factors  suggested  by  the  responsible  construction  agency. _ 
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Figure  A 1.  Instrumentation  for  a  typical  field  site,  plan  view. 


reliable  data  for  correlation  between  field  and  laboratory 
studies.  It  is  recognized  that  other  systems  could  also 
be  feasible.  Forthcoming  technological  developments 
may  provide  less  expensive  and/or  more  reliable  methods 
of  accomplishing  some  observations.  Table  A2  lists 
the  suggested  types  of  instrumentation,  the  vertical 
spacing  of  individual  sensors,  the  frequency  of  obser¬ 
vations  and  the  type  of  data  collection  system.  Table 


A3  lists  the  approximate  cost  of  instrumentation  for 
a  typical  site.  Samples  of  nearly  every  material  in  the 
pavement  profile  must  be  obtained  and  returned  to  the 
coordinating  organization  for  distribution  to  agencies 
that  will  conduct  the  :  boratory  tests.  Similar  soils 
should  be  used  to  replace  those  removed  from  the  pave¬ 
ments  for  the  laboratory  studies. 
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Figure  A  2.  Instrumentation  for  a  typical  field  site,  profile  view. 

Table  A2.  Instrumentation  for  field  sites. 


Item 

Purpose 

Observations 

Thermocouples 

Temperatures 

Every  6  hours 

(automatically) 

Tensiometers 

Stress  in 

Every  6  hours 

soli  water 

(automatically) 

Frost  heave 

Subsurface 

Daily 

disks 

frost  heave 

(automatically) 

Groundwater 

Groundwater 

Daily 

well 

level 

(automatically) 

Surface 

Frost  heave 

3  times  per  week 

elevation 

at  surface 

(manually) 

points 

Benchmark 

Reference 

- 

point 

Instrument 

House 

- 

shelter 

instruments 

_ Spacing _ 

Top  and  bottom  of  pavement,  every  6  in.  (1S.2  cm)  through  frost  zone 
additional  thermocouples  6,  18  and  36  in.  (15.2,  45.7  and  91.4  cm) 
below  design  seasonal  frost  depth. 

Every  6  in.  (15.2  cm)  in  fine  sands  and  fine-grained  soils  in  frost  zone, 
one  6  in.  below  frost  zone  and  one  6  in.  below  water  table. 

Every  6  in.  (15.2  cm)  to  a  depth  of  12  in.  (45.7  cm)  below  subgrade 
surface. 

One  well,  to  a  depth  about  24  in.  (6 1 .0)  below  groundwater  table. 

Starting  from  center  of  test  pit,  lay  out  5X5  ft  (1.52X1.52  m)  grid  for 
20  ft  (6.1  m)  each  longitudinal  direction  and  10  ft  (3.05  m)  each  trans¬ 
verse  direction. 

One  bench  mark  which  will  not  be  moved  by  frost  heave. 

One  digital  data  collection  system  capable  of  scanning  the  desired 
number  of  points  at  the  desired  intervals.  Instrument  shelter  must 
maintain  70°±20°F  (21.1°±1 1.1°C)  and  contain  electrical  power  for 
data  system  and  other  instrumentation  and  heat  and  air  conditioning 
systems. _ 


Table  A3.  Approximate  cost  of  instrumentation  at  a  typical  field  site. 


Thermocouples  -  12  sensors 

Wire  $0.75/ftX  100  ft  ($2.46/mX  30.48  m)  $  75 

Labor  (assembly)  8  hx$10/h  80 

Subtotal 

Tensiometers  -  7  sensors 

Material  SI  600 

Labor  (assemble,  fill  and  calibrate)  40  hx$10/h  400 

Subtotal 

Frost  heave  disks  -  6  sensors 

Material  (6  eaX $20)  $  120 

Labor  (calibration  and  add  extension  wire) 

16  hx$10/h  160 

Read-out  device  2000 

Subtotal 

Groundwater  well  -  1  sensor 

Material  $  10 

Labor  (fabricate)  4  hX$10/h  40 

Subtotal 

Surface  elevation  points  $  10 

Benchmark  100 

Data  system  (30  channels)  7000 

Instrument  shelter  30  ft1  (2.79  m2 )  @  $30/ft’ )  ($322.92/mJ )  900 

Test  pit  (digging  and  refilling)  ( 1 6  hX  $  1 0/h)  160 

Installing  and  connecting  instrumentation  (40hX$10/h)  400 

Drill  rig  and  operators  ($50+8  hx  $  1 0/h)  1 30 

Subtotal 
Total  Cost 


$  155 


2,000 


2,280 


50 


8,700 

$13,275 


Note  that  this  total  does  not  include  indirect  costs,  overhead  costs,  monthly  costs  for  power,  heat 
and  air  conditioning  of  instrument  shelter  or  labor  costs  of  making  observations  manually.  Nor  does 
it  include  costs  for  conducting  laboratory  frost-susceptibility  tests  and  other  laboratory  studies  and 
tests.  If  materials  for  several  test  sites  were  purchased  under  a  single  order,  discounts  on  some  items 
could  be  obtained  for  quantity  lots. 
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APPENDIX  B:  PROPOSED  INVESTIGATION  OF 
THAW  WEAKENING  OF  SUBGRADE  SOIL  AND 
GRANULAR  UNBOUND  BASE  COURSE. 

Submitted  to  the  Federal  Highway  Administration  on 
13  April  1978. 

Most  methods  of  design  of  pavements  whose  sub¬ 
grades  are  affected  by  frost  action  make  use  of  a  single 
index  of  subgrade  supporting  capacity.  Whether  the 
index  be  a  CBR-value,  an  R-value,  a  K-value,  a  soil  sup¬ 
port  value,  a  soil  classification  index,  or  some  other 
term,  it  is  usual  to  assign  for  each  soil  a  single  numerical 
value  of  the  index  that  is  construed  to  be  applicable  in 
design  of  pavements  to  serve  year-round  traffic.  The 
index  in  some  cases  is  selected  to  represent  the  worst 
case,  when  the  subgrade  is  severely  weakened  by  thawing, 
while  in  other  cases  it  does  not  account  for  freezing  and 
thawing  and  is  more  representative  of  the  conditions  in 
summer  or  fall  when  subgrades  have  recovered  from 
the  thaw-weakened  state. 

Use  of  a  single  index  occasions  at  least  two  sources 
of  uncertainty.  On  one  hand  the  design  equations, 
which  currently  are  in  almost  every  case  empirical,  must 
include  among  their  coefficients  some  means  of  ac¬ 
counting  for  the  extreme  seasonal  variation  in  subgrade 
supporting  capacity.  And  on  the  other  hand,  the 
selected  soil  index  and  its  numerical  values  must  be 
able  to  consistently  evaluate  a  great  diversity  of  soils 
having  markedly  different  seasonal  variations  in  prop¬ 
erties.  In  view  of  these  uncertainties,  it  seems  hardly 
surprising  that  different  pavements  constructed  in  frost 
areas  vary  from  conditions  with  no  evidence  of  sub¬ 
grade-related  distress  (indicating  good  design  or,  one 
might  suppose,  even  overdesign),  to  those  with  exten¬ 
sive  evidence  of  subgrade-related  distress  (indicating 
underdesign). 

As  mechanistic  design  methods  (those  based  on 
calculated  stresses  and  strains)  become  more  widely 
adopted,  it  will  be  possible  to  dispel  the  first  uncertainty 
by  incorporating  into  the  method  a  damage  accumulator. 
The  pavement  damage  that  will  occur  in  each  season 
can  be  calculated  and  summed  to  represent  the  total 
yearly  damage  occasioned  by  traffic.  A  trial  pavement 
design  would  then  be  adopted  for  construction  if  the 
total  damage  accumulation  during  the  selected  number 
of  years  of  economic  life  were  within  limits  determined 
in  relation  to  desired  pavement  serviceability.  Effective 
application  of  a  cumulative  damage  approach,  however, 
requires  that  the  actual  properties  of  subgrade  soils  in 
each  season  be  determined;  accomplishment  of  this 
task  will  also  dispel  the  second  uncertainty  mentioned 
above. 

Design  procedures  are  available  (Bergan  and  Moni- 
smith  1972,  Barker  and  Brabston  1975)  for  analysis  of 


damage  accumulation  in  pavements  subjected  to  freezing 
and  thawing.  Most  methods  require  a  prediction  of  the 
seasonal  values  of  the  resilient  modulus  of  each  type  of 
soil,  base  course,  or  other  material  in  the  pavement 
profile.  While  certain  information  is  available  (Johnson 
et  al.  1975)  regarding  effects  of  frost  action  on  the 
resilient  modulus  of  soils,  further  research  is  needed  to 
develop  laboratory  test  techniques  for  determining 
seasonal  values  of  the  resilient  modulus  for  design  of 
new  pavements  and  rehabilitation  of  existing  pavements. 
The  resilient  properties  need  to  be  determined  and  ex¬ 
pressed  in  relation  to  pertinent  environmental  parameters 
or  to  soil  properties  that  are  sensitive  to  changes  in  such 
environmental  conditions.  It  is  expected  that  the  state 
of  stress  in  the  soil  moisture  will  prove  to  be  a  con¬ 
venient  parameter,  as  previous  research  (Fredlund  et  al. 
1975)  has  shown  that  one  of  the  most  significant  vari¬ 
ables  affecting  the  resilient  modulus  is  the  moisture 
tension.  Moisture  tension  has  the  further  powerful 
advantage  that  with  the  recent  development  of  frost- 
resistant  tensiometers  (McKim  et  al.  1976)  it  can  be 
monitored  continuously  in  situ. 

There  is  currently  a  critical  need  for  a  broad  data  base 
of  actual  values  of  the  resilient  modulus  for  a  wide 
variety  of  types  of  subgrade  soils  and  base  course  materials, 
measured  throughout  the  year  and  including  periods  of 
freezing,  thawing,  and  recovery.  Lack  of  information  on 
seasonal  changes  in  the  resilient  modulus  constitutes  an 
impediment  to  the  implementation  of  mechanistic  design 
methods.  Such  a  data  base  would  be  of  immediate  and 
direct  use  in  pavement  design  and,  with  the  corresponding 
soil  and  environmental  data,  could  be  also  used  to  develop 
a  general  model  of  seasonal  changes  in  resilient  moduli 
of  subgrade  and  base  course  soils. 

The  planned  field  investigations  for  validation  of 
soil  frost-susceptiblity  tests  provide  an  excellent  op- 
protunity  for  collection  of  resilient  modulus  data  on  a 
number  of  subgrade  soils  and/or  base  courses  of  various 
types  and  properties.  We  propose  to  conduct  both 
field  and  laboratory  tests  for  that  purpose. 

A  repeated-load  plate  bearing  test  is  recommended 
for  the  field  test.  In  this  test  repeated  loads  are  applied 
to  the  test  pavement  through  a  circular  plate  or  other 
loading  medium  and  the  recoverable,  or  resilient,  de¬ 
formation  is  recorded  at  the  plate  and  on  the  adjacent 
pavement  at  set  radial  distances  from  the  plate.  The 
magnitude,  frequency  and  wave-form  of  the  cyclic  loads 
are  established  to  typify  selected  critical  traffic  con¬ 
ditions  on  highways  or  airfields.  Resilient  moduli  of 
certain  materials  in  the  layered  pavement  system  are 
determined  by  separate  laboratory  tests  and  are  input, 
together  with  data  from  the  plate-bearing  tests,  to  a 
selected  pavement  response  model  to  determine  the 
resilient  modulus  of  the  subgrade  or  base  course  layer  in 
question. 
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We  have  a  project  currently  in  progress  at  CRREL 
to  obtain  data  in  the  manner  outlined  above.  To  date 
we  have  tested  one  silt  and  one  clay  subgrade  soil 
Johnson  et  al.  (1976)  and  have  recently  completed  the 
first  freeze-thaw  cycle  in  tests  on  a  silty  sand.  We  have 
found  it  convenient  to  employ  "full-depth,"  or  all- 
bituminous-concrete  (ABC)  pavements,  to  minimize 
the  number  of  layers  in  the  system.  Thus  far,  we  have 
been  able  to  use  a  5-layer  solution  (CHEV-5L),  although 
in  some  cases  the  system  has  required  use  of  a  larger 
modified  Chevron  program,  which  admits  up  to  15 
layers.  The  BISTRO*  program  also  has  been  used,  and 
both  the  BISAR*  program  and  a  finite  element  model 
also  are  available  and  could  be  used  if  convenient. 

The  pavements  selected  for  repeated-load  plate 
bearing  tests  would  be  the  same  ones  used  for  the  frost- 
susceptibility  validation  testing  mentioned  above,  which 
forms  part  of  the  same  research  project  directed  toward 
pavements  in  frost  areas.  For  both  types  of  tests  it  would 
be  desirable  to  select  test  pavements  at  field  sites  having 
a  wide  range  in  subgrade  soil  types  and  moisture  con¬ 
ditions.  Also,  for  better  effectiveness  in  achieving  the 
objectives  of  both  types  of  tests,  it  is  desirable  that  the 
sites  have  simple  pavement  sections,  since  the  com¬ 
plexity  of  multi-layered  sections  makes  the  acquisition 
of  data  on  individual  soils  more  difficult.  In  particular, 
it  would  be  preferable  that  several  sites  have  ABC  pave¬ 
ments  or  pavements  with  only  one  type  of  base  course 
material. 

The  instrumentation  planned  to  be  installed  at  the 
test  sites  for  frost-susceptibility  validation  tests  would 
be  sufficient  also  for  thaw-weakening.  At  each  of  the 
pavement  sections  selected  for  resilient  modulus  testing 
we  propose  to  perform  6  to  10  loading  tests  during  one 
calendar  year  for  determination  of  the  resilient  modulus 
of  the  subgrade  soil  or  base  material. 

The  first  loading  test  would  be  conducted  when  the 
sections  were  deeply  frozen,  followed  by  tests  when 
thawing  advanced  into  the  layer  in  question  about  5, 

10,  20,  and  40  cm.  Further  tests  would  be  made  when 
the  soil  profile  had  just  reached  a  fully  thawed  state, 
and  then  after  elapsed  times  of  about  1 5,  30,  60  and 
100-200  days.  Depending  on  the  conditions,  certain 
of  the  tests  might  be  omitted.  When  each  test  was  per¬ 
formed,  temperatures  and  moisture  tension  would  be 
recorded  by  means  of  the  installed  thermocouples  and 
tensiometers,  and  in  some  cases,  a  shallow  boring  would 
be  made  nearby  to  obtain  moisture  contents  in  the  pave¬ 
ment  section  and  subgrade.  Prior  to  freezing  and  again 
when  the  sections  were  deeply  frozen,  special  borings 
would  be  taken  to  obtain  undisturbed  cores  of  the 
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base  and  subgrade  materials  and  the  bituminous  con¬ 
crete  pavement  for  laboratory  testing.  The  data  from 
each  plate  bearing  test  would  be  analyzed,  by  means  of 
the  selected  pavement  response  model,  to  calculate  the 
resilient  modulus  of  the  base  course  or  subgrade  soil 
in  question.  The  entire  test  series  at  each  site,  accordingly, 
would  show  the  variation  in  modulus  of  that  soil  through¬ 
out  a  complete  annual  cycle.  Of  special  interest  are  not 
only  the  absolute  values  of  resilient  moduli  and  their 
relationship  with  moisture  content  and  moisture  tension, 
but  also  the  duration  of  the  thaw-weakened  condition 
and  the  rate  of  subsequent  recovery. 

It  would  also  be  useful  to  perform  conventional 
static  load  testing  at  some  of  the  test  sites  to  determine 
the  CBR  and  subgrade  modulus  (k,  psi/in.)  of  the  various 
materials.  As  these  are  destructive  tests  requiring  excava¬ 
tion  of  pits  to  perform  the  tests  on  the  surface  of  each 
layer,  the  tests  probably  would  be  feasible  at  only  a 
few  sites.  Wherever  such  tests  could  be  performed,  it 
would  be  useful  to  compare  the  data  with  the  resilient 
moduli  determined  by  field  and  laboratory  testing. 

We  propose  to  make  use  of  the  CRREL  repeated- 
load  plate  bearing  test  device  (a  large  semitrailer-mounted 
rig)  at  any  sites  located  within  a  reasonable  distance  from 
Flanover,  New  Hampshire— possibly  a  radial  distance  of 
about  150  miles.  At  other  sites  it  is  hoped  the  FHWA 
test  van  might  be  made  available,  and  at  more  distant 
sites  other  tests  will  be  used,  depending  on  the  equip¬ 
ment  available  in  the  locality.  One  possibility  is  Dyna- 
flect  equipment  or  other  similar  devices,  and  the  Benkel- 
man  beam  could  also  be  used.  These  devices  are  far  from 
ideal,  however,  in  the  first  case  because  of  the  low  stress 
levels  and  in  the  second  because  of  the  low-speed,  non- 
repetitive  loading.  Consequently,  if  no  better  alternatives 
are  available,  we  would  expect  that  at  some  sites  it  would 
be  considered  acceptable  to  forgo  field  loading  tests  and 
to  determine  the  resilient  properties  of  the  materials 
exclusively  by  means  of  laboratory  tests.  In  these  cases 
the  seasonal  change  of  moisture  content  and  moisture 
tension  in  the  various  layers  of  the  pavement  test  sec¬ 
tion  would  be  monitored  to  make  it  possible  to  relate 
the  laboratory  test  moduli  to  specific  times  during  the 
recovery  process. 

To  obtain  undisturbed  core  samples  of  unfrozen  and 
frozen  subgrade  and  base  course  materials  it  would  be 
convenient  to  use  CRREL  drilling  equipment,  except 
at  distant  sites  where  equipment  available  in  the  locality 
would  be  used.  In  some  pavement  sections,  we  expect 
that  it  would  prove  impossible  to  obtain  intact  core 
samples  of  certain  of  the  coarser-grained  materials  that 
may  be  poorly  bonded  even  in  the  frozen  state,  due  to 
a  low  ice  content  or  to  temperatures  only  slightly  below 
freezing.  In  these  cases  loose  samples  would  be  taken, 
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and  remolded  specimens  would  be  prepared  and  frozen 
in  the  laboratory.  If  no  better  alternative  is  found,  it 
may  be  decided  in  some  cases  to  forgo  laboratory  re¬ 
silient  modulus  testing  and  rely  exclusively  on  field 
plate  bearing  tests  for  determination  of  seasonal  variation 
in  resilient  moduli  of  the  various  materials. 

The  associated  laboratory  testing  program  would 
utilize  a  repeated  load  triaxial  system.  Subgrade  and 
base  course  materials  obtained  from  each  test  site  in 
the  frozen  state  would  be  subjected  to  controlled  re¬ 
peated  triaxial  loading  in  the  frozen  and  thawed  condi¬ 
tions  to  obtain  resilient  moduli  and  Poisson's  ratios.  In 


addition,  undisturbed  or  remolded  unfrozen  samples  i 

would  be  tested  to  evaluate  the  range  of  seasonal  con-  j 

ditions.  Thawed  samples  would  be  tested  over  a  range  j 

of  moisture  contents  in  order  to  simulate  the  recovery  j 

period.  Attempts  would  be  made  to  measure  the  pre-  I 

vailing  stress  state  in  the  pore  water  before  each  series  \ 

of  repeated  loads  is  applied.  j 


Because  the  resilient  modulus  of  an  asphalt  concrete  j 

pavement  changes  with  temperature  and  moisture  con-  I 

tent,  laboratory  repeated-load  tests  would  also  be  con-  j 

ducted  at  various  environmental  conditions  on  pave-  i 

ment  cores  obtained  from  the  field  sites.  The  applied 
load  duration  and  frequency  would  duplicate  that  ap¬ 
plied  in  the  repeated-load  plate  bearing  test.  The  tests  j 

would  be  conducted  over  a  range  of  deviator  and  con-  ■ 

fining  stress  amplitudes  to  cover  the  range  of  stress  i 

levels  anticipated  at  the  test  sites.  The  test  results  would  j 

be  evaluated  by  statistical  techniques,  the  product  being  I 

the  resilient  modulus  and  Poisson’s  ratio  as  a  function  of 
the  important  variables.  We  expect  that  moisture  ten¬ 
sion  would  be  a  principal  parameter  in  samples  in  which 
it  can  be  successfully  measured  and  continuously  moni¬ 
tored.  Other  variables  to  be  evaluated  include  the  initial  j 

state  (frozen,  thawed,  or  never  frozen),  water  content,  j 

density  and  the  deviator  and  confining  stress  levels.  The 

analysis  selects  the  significant  variables  and  provides  co-  S 

efficients  necessary  to  calculate  the  seasonal  values  of  the 

resilient  modulus  or  Poisson's  ratio.  The  results,  thus  I 

determined,  would  be  compared  with  the  results  obtained  j 

from  the  analysis  of  the  repeated-load  plate  bearing 

tests  on  the  field  sections.  j 

An  attempt  would  also  be  made  to  develop  a  general  j 

model  for  the  resilient  modulus  of  the  various  elements  j 

of  a  pavement  system.  It  is  anticipated  that  this  would 

be  an  empirical  model  based  on  observed  relationships  ] 

between  the  resilient  modulus  and  parameters  such  as 

moisture  tension,  moisture  content,  density,  tern-  j 

perature,  permeability,  time,  confining  stress,  Poisson’s  : 

ratio,  and  indices  of  material  type.  The  success  of  the  j 

frost  heave  model  would  be  dependent  on  the  capability  | 

of  this  model  to  predict  heave,  and  thus  changes  in  1 

moisture  content  and  density  under  field  conditions.  ] 
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APPENDIX  C:  DERIVATION  OF  FINITE  ELEMENT  SYSTEM  MATRICES 


In  this  appendix,  the  solution  of  eq  41  is  carried  out  term  by  term  for  the  mth  interior  element. 
In  this  derivation  x  will  be  considered  the  local  coordinate. 

Let  [s]  be  a  3  by  3  square  matrix  and  be  defined  by  the  first  two  terms  of  eq  43  as  follows  and 
let  <N>  be  defined  as  in  the  main  text.  In  this  derivation,  eq  38  is  expressed  as: 
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Let  [p]  be  a  3  by  3  square  matrix,  <N>  be  defined  as  previously  and  let  [p]  equal  the  third  term 
in  eq  43  as  follows: 
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Let {f}  be  a  column  matrix,  and  let  <N>  be  defined  as  previously  given,  and  let  {fl  equal  the  last 
term  in  eq  43  modified  as  follows: 
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Finally,  assembling  all  the  above  into  standard  matrix  notation 
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These  equations  are  the  element  system  equation  in 

[si  (u)  +  (p  1  (ul  -  {£} 

which  is  the  finite  element  solution  of  eq  43  modified  for  the  spatial  domain  only.  The  temporal 
solution  is  by  the  Crank-Nicholson  method  and  was  discussed  in  the  text.  The  following  section  pre¬ 
sents  an  assemblage  of  the  M  element  matrices  into  the  solution  scheme  for  the  entire  solution 
domain. 

Assembly  of  system  matrices 

The  global  system  matrices  are  obtained  by  appropriately  assembling  the  element  solutions  ob¬ 
tained  before.  After  assemblage  the  result  may  be  generalized  as  follows: 

(s  I  (u)  +  [?  1  {u;  -  {Fl 

where  [S]  is  a  square  banded  nonsymmetrical  matrix  because  of  the  k2  term  in  [s],  [P]  is  a  square 
banded  symmetric  matrix  and  {F}  is  a  column  matrix.  The  matrices  are  functions  of  the  k:  para¬ 
meters  and  geometrical  discretion  of  the  system  and  are  known.  The  unknown  matrices  are  u  and 
{ujwhich  are  the  nodal  values  of  the  state  variable  for  each  node  of  the  problem.  Recall  that  the 
dot  (in  (u })  refers  to  the  time  derivative. 

The  system  matrices  are  as  follows  before  considering  the  boundary  conditions: 
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[Pj  is  identically  formed,  and 


The  superscript  refers  to  the  element  number  and  the  subscript  refers  to  the  values  position  in  the 
element  system  matrices. 

Considering  specified  boundary  conditions  at  the  ends  of  the  solution  domain,  ux=0  =  uu  and 
ux=L  =  UL> the  above  matrices  became 


0 


0 


If  a  boundary  is  assumed  to  be  a  natural  boundary,  i.e.  3u/dx  =  0,  modification  of  the  system 
matrices  is  not  required.  For  instance,  assume  that  at  x  =  0  there  is  a  natural  boundary  condition, 
then  the  first  part  of  the  system  matrices  would  be  the  same  as  the  first  unmodified  matrices  above, 
A  second  order  polynomial  (quadratic)  shape  function  was  used  to  develop  the  finite-element 
analog.  Other  shape  functions  are  also  possible.  Previous  work,  Guymon  (1976),  contained  a 
derivation  of  the  element  systems  matrices  for  a  linear  shape  function.  The  use  of  a  second  order 
shape  function  requires  more  core  storage  for  computer  solution;  however,  accuracy  has  been 
improved. 

Computer  storage  of  arrays 

A  substantial  benefit  of  the  finite-element  method  is  that  the  system  matrices  are  '‘naturally” 
banded.  Consequently,  computer  storage  arrays  can  be  substantially  reduced  if  only  the  nonzero 
information  in  the  arrays  is  stored.  This  is  accomplished  as  is  illustrated  in  the  following  figure: 
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Thus,  if  the  number  of  nodes  in  a  problem  is  N,  the  required  storage  for  the  [Sj  array  is  N  by  5.  If 
the  entire  array  were  stored,  the  storage  required  would  be  N  by  N.  Algorithms  are  specially  coded 
to  accommodate  the  banded  storage  concept  which  also  saves  substantial  computer  execution  time 
since  operations  with  zeros  are  not  carried  out. 
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A  facsimile  catalog  card  in  Library  of  Congress  MARC 
format  is  reproduced  below. 
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